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All galaxies once passed through a hyper-luminous quasar phase powered by 
accretion onto a supermassive black hole. But because these episodes are brief, 
quasars are rare objects typically separated by cosmological distances. In a 
survey for Lyman-o emission at redshift z & 2, we discovered a physical as¬ 
sociation of four quasars embedded in a giant nebula. Located within a sub¬ 
stantial overdensity of galaxies, this system is likely the progenitor of a massive 
galaxy cluster. The chance probability of finding a quadruple-quasar is esti¬ 
mated to be ~ Hr 7 , implying a physical connection between Lyman-a nebulae 
and the locations of rare protoclusters. Our findings imply that the most mas¬ 
sive structures in the distant Universe have a tremendous supply (~ 10 n solar 
masses) of cool dense (volume density ~ 1 cm -3 ) gas, in conflict with current 
cosmological simulations. 

Cosmologists do not fully understand the origin of supermassive black holes (SMBHs) at the 
centers of galaxies, and how they relate to the evolution of the underlying dark matter, which 
forms the backbone for structure in the Universe. In the current paradigm, SMBHs grew in 
every massive galaxy during a luminous quasar phase, making distant quasars the progenitors 
of the dormant SMBHs found at the centers of nearby galaxies. Tight correlations between the 
masses of these local SMBHs and both their host galaxy ( 1 ) and dark matter halo masses (2) 
support this picture, further suggesting that the most luminous quasars at high-redshift should 
reside in the most massive galaxies. It has also been proposed that quasar activity is triggered 
by the frequent mergers that are a generic consequence of hierarchical structure formation (3,4). 
Indeed, an excess in the number of small-separation binary quasars (5, 6), as well as the mere 
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existence of a handful of quasar triples (7,8), support this hypothesis. If quasars are triggered by 
mergers, then they should preferentially occur in rare peaks of the density field, where massive 
galaxies are abundant and the frequency of mergers is highest (9). 

Following these arguments, one might expect that at the peak of their activity z ~ 2 — 3, 
quasars should act as signposts for protoclusters, the progenitors of local galaxy clusters and 
the most massive structures at that epoch. However, quasar clustering measurements (10,11) 
indicate that quasar environments at z ~ 2 — 3 are not extreme: these quasars are hosted by 
dark matter halos with masses M halo ~ 10 12 5 M 0 (where M e is the mass of the Sun), too 
small to be the progenitors of local clusters (12). But the relationship between quasar activity 
and protoclusters remains unclear, owing to the extreme challenge of identifying the latter at 
high-redshift. Indeed, the total comoving volume of even the largest surveys for distant galaxies 
at z ~ 2 — 3 is only ~ 10' Mpc 3 , which would barely contain a single rich cluster locally. 

Protoclusters have been discovered around a rare-population of active galactic nuclei (AGN) 
powering large-scale radio jets, known as high-redshift radio galaxies (HzRGs) (13). The 
HzRGs routinely exhibit giant ~ 100 kpc nebulae of luminous Lyo emission L Ly a ~ 10 44 ergs 4 . 
Nebulae of comparable size and luminosity have also been observed in a distinct population of 
objects known as ‘Lycc blobs’ (LABs) (14,15). The LABs are also frequently associated with 
AGN activity (16-18), although lacking powerful radio jets, and appear to reside in overdense, 
protocluster environments (15,19, 20). Thus among the handful of protoclusters (21) known, 
most appear to share two common characteristics: the presence of an active SMBH and a giant 
Lya nebula. 

We have recently completed a spectroscopic search for extended Lycc emission around a 
sample of 29 luminous quasars at z ~ 2, each the foreground (f/g) member of a projected 
quasar pair (22). Analysis of spectra from the background (b/g) members in such pairs reveals 
that quasars exhibit frequent absorption from a cool, metal-enriched, and often optically thick 
medium on scales of tens to hundreds kpc (22-28). The UV radiation emitted by the luminous 
f/g quasar can, like a flashlight, illuminate this nearby neutral hydrogen, and power a large-scale 
Lya-cmission nebula, motivating our search. By construction, our survey selects for exactly 
the two criteria which seem to strongly correlate with protoclusters: an active SMBH and the 
presence of a large-scale Lyo emission nebula. 

Of the 29 quasars surveyed, only SDSSJ0841+3921 exhibited extended large-scale (> 50 kpc) 
Lycc emission above our characteristic sensitivity limit of 6 x 10~ 18 erg s -1 cm -2 arcsec -2 (2a). 
We designed a custom narrow-band filter tuned to the wavelength of Lyo at the f/g quasar red- 
shift z = 2.0412 (A C e n ter = 3700A, FWHM A = 33A), and imaged the field with the Keck/LRIS 
imaging spectrometer for 3hrs on UT 12 November 2012. The combined and processed images 
reveal Lya emission from a giant filamentary structure centered on the f/g quasar and extending 
continuously towards the b/g quasar (see Fig. 1). This nebulosity has an end-to-end size of 37" 
corresponding to 310 kpc and a total line luminosity ^Lya = 2.1 x 10 44 ergs _1 , making it one 
of the largest and brightest nebula ever discovered. 

The giant nebula is only one of the exceptional properties of SDSSJ 08414-3921. Our 
images reveal three relatively compact candidate Lyct emitting sources with faint continuum 


2 



magnitudes V ~ 23 — 24, embedded in the Lyct filament and roughly aligned with its major 
axis. Follow-up spectroscopy reveals that the sources labeled AGN1, AGN2, and AGN3 are 
three AGN at the same redshift as the f/g quasar (see right panel of Fig. land (29)), making this 
system the only quadruple AGN known. Adopting recent measurements of small-scale quasar 
clustering (30), we estimate that the probability of finding three AGN around a quasar with 
such small separations is ~ 10~ 7 (31). Why then did we discover this rare coincidence of AGN 
in a survey of just 29 quasars? Did the giant nebula mark the location of a protocluster with 
dramatically enhanced AGN activity? 

To test this hypothesis, we constructed a catalog of Lyo;-emitting galaxies (LAEs), and 
computed the cumulative overdensity profile of LAEs around SDSSJ0841+3921, relative to 
the background number expected based on the LAE luminosity function (32) (see Fig. 2). To 
perform a quantitative comparison to other giant Lya nebulae, many of which are known to 
coincide with protoclusters, we measured the giant nebulae-LAE cross-correlation function for 
a sample of eight systems - six HzRGs and two LABs - for which published data was available 
in the literature (33). In Fig. 2, we compare the overdensity profile around SDSSJ0841+3921 to 
this giant nebulae-LAE correlation function. On average, the environment of HzRGs and LABs 
hosting giant Lya nebulae (red line) is much richer than that of radio-quiet quasars (10) (blue 
line), confirming that they indeed reside in protoclusters. Furthermore, the clustering of LAEs 
around SDSSJ0841+3921 has a steeper overdensity profile, and exceeds the average protoclus¬ 
ter by a factor of > 20 for R < 200 kpc and by ~ 3 on scales of R ~ 1 Mpc. In addition to 
the overdensity of four AGN, the high number of LAEs surrounding SDSSJ0841+3921 make 
it one of the most overdense protoclusters known at z ~ 2 — 3. 

The combined presence of several bright AGN, the Lyo: emission nebula, and the b/g absorp¬ 
tion spectrum, provide an unprecedented opportunity to study the morphology and kinematics 
of the protocluster via multiple tracers, and we find evidence for extreme motions (34). Specif¬ 
ically, AGN 1 is offset from the precisely determined systemic redshift (35) of the f/g quasar by 
+1300 ± 400 km s’ 1 . This large velocity offset cannot be explained by Hubble expansion - the 
miniscule probability of finding a quadruple quasar in the absence of clustering P rsj 10~ 13 (31) 
and the physical association between the AGN and giant nebula demand that the four AGN 
reside in a real collapsed structure - and thus provides an unambiguous evidence for extreme 
gravitational motions. In addition, our slit spectra of the giant Lya nebula reveal extreme kine¬ 
matics of diffuse gas (Fig. 3), extending over a velocity range of —800 to +2500 km s -1 from 
systemic. Furthermore, there is no evidence for double-peaked velocity profiles characteristic 
of resonantly-trapped Lya, which could generate large velocity widths in the absence of corre¬ 
spondingly large gas motions. Absorption line kinematics of the metal-enriched gas, measured 
from the b/g quasar spectrum at impact parameter of R± = 176 kpc (Fig. 4), show strong ab¬ 
sorption at « +650 km s^ 1 with a significant tail to velocities as large as ~ 1000 km s -1 . It is 
of course possible that the extreme gas kinematics, traced by Lyo emission and metal-line ab¬ 
sorption, are not gravitational but rather arise from violent large-scale outflows powered by the 
multiple AGN. While we cannot completely rule out this possibility, the large velocity offset of 
+1300 km s’ 1 between the f/g quasar and the emission redshift of AGN1 can only result from 
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gravity. 

One can only speculate about the origin of the dramatic enhancement of AGN in the 
SDSSJ0841+3921 protocluster. Perhaps the duty cycle for AGN activity is much longer in pro¬ 
toclusters, because of frequent dissipative interactions (5,6), or an abundant supply of cold gas. 
A much larger number of massive galaxies could also be the culprit, as AGN are known to trace 
massive halos at z ~ 2. Although SDSSJ0841+3921 is the only example of a quadruple AGN 
with such small separations, previously studied protoclusters around HzRGs and LABs, also 
occasionally harbor multiple AGN (13,36). Regardless, our discovery of a quadruple AGN and 
protocluster from a sample of only 29 quasars suggests a link between giant Lya nebulae, AGN 
activity, and protoclusters - similar to past work on HzRGs and LABs - with the exception that 
SDSSJ 0841+3921 was selected from a sample of normal radio-quiet quasars. From our survey 
and other work (37), we estimate that ~ 10% of quasars exhibit comparable giant Lya nebulae. 
Although clustering measurements imply that the majority of z ~ 2 quasars reside in moderate 
overdensities (10-12), we speculate that this same 10% trace much more massive protoclus¬ 
ters. SDSSJ0841+3921 clearly supports this hypothesis, as does another quasar-protocluster 
association (10,38), around which a giant Lya nebula was recently discovered (39,40). 

Given our current theoretical picture of galaxy formation in massive halos, an association 
between giant Lya nebulae and protoclusters is completely unexpected. The large Lya lumi¬ 
nosities of these nebulae imply a substantial mass ( rs./ 10 11 M 0 ) of cool (T ~ 10 4 K) gas (41), 
whereas cosmological hydrodynamical simulations indicate that already by z ~ 2 — 3, baryons 
in the massive progenitors (M ha i 0 > 10 13 M 0 ) of present-day clusters are dominated by a hot 
shock-heated plasma T ~ 10' K (42,43). These hot halos are believed to evolve into the X-ray 
emitting intra-cluster medium observed in clusters, for which absorption-line studies indicate 
negligible < 1% cool gas fractions (44). Clues about the nature of this apparent discrepancy 
come from our absorption line studies of the massive ~ 10 12 ' 5 M e halos hosting z ~ 2 — 3 
quasars. This work reveals substantial reservoirs of cool gas > 10 10 M e (22-28), manifest as 
a high covering factor ~ 50% of optically thick absorption, several times larger than predicted 
by hydrodynamical simulations (42,43). This conflict most likely indicates that current simula¬ 
tions fail to capture essential aspects of the hydrodynamics in massive halos at z ~ 2 (27,42), 
perhaps failing to resolve the formation of clumpy structure in cool gas (41). 

If illuminated by the quasar, these large cool gas reservoirs in the quasar circumgalactic 
medium (CGM) will emit fluorescent Lya photons, and we argue that this effect powers the neb¬ 
ula in SDSSJ0841+3921 (45). But according to this logic, nearly every quasar in the Universe 
should be surrounded by a giant Lya nebulae with size comparable to its CGM (~ 200 kpc). 
Why then are these giant nebulae not routinely observed? 

This apparent contradiction can be resolved as follows. If this cool CGM gas is illumi¬ 
nated and highly ionized, it will fluoresce in the Lya line with a surface brightness scaling 
as SB LyQ oc A H nn, where Ah is the column density of cool gas clouds which populate the 
quasar halo, and % is the number density of hydrogen atoms inside these clouds. Note the total 
cool gas mass of the halos scales as M cool oc f? 2 A H , where R is the radius of the halo (45). 
Given our best estimate for the properties of the CGM around typical quasars (% ~ 0.01 cm -3 
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and A'n ~ 10 2 ° cm 2 or M coo \ ~ 10 11 M 0 ) (22, 26, 27), we expect these nebulae to be ex¬ 
tremely faint SB Lya ~ 10~ 19 ergs -1 cm -2 arcsec -2 , and beyond the sensitivity of current in¬ 
struments (22). One comes to a similar conclusion based on a full radiative transfer calculation 
through a simulated dark-matter halo with mass M hak) « 10 12,5 M 0 (41). Thus the factor of 
~ 100 times larger surface brightness observed in the SDSSJ0841+3921 and other protocluster 
nebulae, arises from either a higher %, iV H (and hence higher M coo 0, or both. The cool gas 
properties required to produce the SDSSJ0841+3921 nebula can be directly compared to those 
deduced from an absorption line analysis of the b/g quasar spectrum (46). 

The b/g quasar sightline pierces through SDSSJ0841+3921 at an impact parameter of R± = 
176 kpc, giving rise to the absorption spectrum shown in Fig. 4. Photoionization modeling of 
these data constrains the total hydrogen column density to be \og l0 Nn = 20.4 ± 0.4 (45), 
implying a substantial mass of cool gas 1.0 x 10 11 M 0 < M cooi < 6.5 x 10 11 M 0 within 
r = 250 kpc. Assuming that the Ly a emitting gas has the same column density as the gas 
absorbing the b/g sightline, reproducing the large fluorescent Lya surface brightness requires 
that this gas be distributed in compact r c i ou d ~ 40 pc clouds at densities characteristic of the 
interstellar medium n\\ ~ 2 cm' 3 , but on ~ 100 kpc scales in the protocluster. 

Clues to the origin of these dense clumps of cool gas comes from their high enrichment level, 
which we have determined from our absorption line analysis (46) to be greater than l/10th of 
the solar value. At first glance, this suggests that strong tidal interactions due to merger activity 
or outflows due to powerful AGN feedback are responsible for dispersing dense cool gas in the 
protocluster. However, the large cool gas mass ~ 10 11 M 0 and high velocities ~ 1000 km s -1 , 
imply an extremely large kinetic luminosity L win( j ~ 10 44 6 for an AGN powered wind, making 
the feedback scenario implausible (25). An even more compelling argument against a merger or 
feedback origin comes from the extremely small cloud sizes r cloud ~ 40 pc implied by our mea¬ 
surements. Such small clouds moving supersonically ~ 1000 km s _1 through the hot T ~ 10' K 
shock-heated plasma predicted to permeate the protocluster, will be disrupted by hydrodynamic 
instabilities in ~ 5 x 10 6 yr, and can thus only be transported ~ 5 kpc (47). These short dis¬ 
ruption timescales instead favor a scenario where cool dense clouds are formed in situ, perhaps 
via cooling and fragmentation instabilities, but are short-lived. The higher gas densities might 
naturally arise if hot plasma in the incipient intra-cluster medium pressure confines the clouds, 
compressing them to high densities (48, 49). Emission line nebulae from cool dense gas has 
also been observed at the centers of present-day cooling flow clusters (50, 51), albeit on much 
smaller scales < 50 kpc. The giant Lya nebulae in z ~ 2 — 3 protoclusters might be manifesta¬ 
tions of the same phenomenon, but with much larger sizes and luminosities, reflecting different 
physical conditions at high-redshift. 

The large reservoir of cool dense gas in the protocluster SDSSJ0841+3921, as well as 
those implied by the giant nebulae in other protoclusters, appear to be at odds with our current 
theoretical picture of how clusters form. This is likely symptomatic of the same problem of too 
much cool gas in massive halos already highlighted for the quasar CGM (27,42,43). Progress 
will require more cosmological simulations of massive halos M > 10 13 M 0 at z ~ 2, as well 
as idealized higher-resolution studies. In parallel, a survey for extended Lya emission around 
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~ 100 quasars would uncover a sample of ~ 10 giant Lya nebulae likely coincident with 
protoclusters, possibly also hosting multiple AGN, and enabling continued exploration of the 
relationship between AGN, cool gas, and cluster progenitors. 
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Fig. 1: Narrow and broad band images of the field surrounding SDSSJ0841+3921 Left: 

Continuum-subtracted, narrow-band image of the field around f/g quasar. The color map and the contours 
indicate the Lya surface brightness (upper color bar) and the signal-to-noise ratio per arcsec 2 aperture 
(lower color bar), respectively. This image reveals a giant Lya nebula on the northern side of the f/g 
quasar and several compact, bright Lya emitters in addition to the f/g quasar. Three of these have been 
spectroscopically confirmed as active galactic nuclei (AGN) at the same redshift. Right: Corresponding 
V -band continuum image of the field presented at left with the locations of the four AGN marked. The 
AGN are roughly oriented along a line coincident with the projected orientation of the Lya nebula. We 
also mark the position of the b/g quasar, which is not physically associated with the quadruple AGN 
system, but whose absorption spectrum probes the gaseous environment of the f/g quadruple AGN and 
protocluster (see Fig. 4). 
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Fig. 2: Characterization of the protocluster environment around SDSSJ0841+3921 . The 

data points indicate the cumulative overdensity profile of LAEs <5(< R) as a function of impact parameter 
R from the f/g quasar in SDSSJ0841+3921, with Poisson error bars. The red curve shows the predicted 
overdensity profile, based on our measurement of the giant nebulae-LAE cross-correlation function de¬ 
termined from a sample of eight systems - six HzRGs and two LABs - for which published data was 
available in the literature. Assuming a power-law form for the cross-correlation £ cross = (r/?’o) -7 , we 
measure the correlation length ro = 29.3 ± 4.9 h ~ 1 Mpc, for a fixed value of 7 = 1.5. The gray shaded 
region indicates the 1 a error on our measurements based on a bootstrap analysis, where both ro and 
7 are allowed to vary. The solid blue line indicates the overdensity of Lyman break galaxies (LBGs) 
around radio-quiet quasars based on recent measurements (10), with the dotted blue lines the ler error 
on this measurement. On average, the environment of HzRGs and LABs hosting giant Lya nebulae is 
much richer than that of radio-quiet quasars (10), confirming that they indeed reside in protoclusters. 
SDSSJ0841+3921 exhibits a dramatic excess of LAEs compared to the expected overdensities around 
radio-quiet quasars (blue curve). Its overdensity even exceeds the average protocluster (red curve) by a 
factor of > 20 for R < 200 kpc decreasing to an excess of ~ 3 on scales of II ~ 1 Mpc, and exhibits a 
much steeper profile. 
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Fig. 3: Ly a spectroscopy of the giant nebula and its associated AGN. Upper Left: The 

spectroscopic slit locations (white rectangles) for three different slit orientations are overlayed on the 
narrow band image of the giant nebula. The locations of the f/g quasar (brightest source), b/g quasar, 
AGN1, and AGN2 are also indicated. Two-dimensional spectra for Slitl (Upper Right), Slit2 (Lower 
Left), and Slit3 (Lower Right) arc shown in the accompanying panels. In the upper right and lower 
left panels, spatial coordinates refer to the relative offset along the slit with respect to the f/g quasar. 
Spectra of AGN1 are present both in Slit 1 (upper right) and Slit 3 (lower right) at spatial offsets 75kpc 
and 25kpc, respectively, while the Lya spectrum of AGN2 is located at a spatial offset — 60kpc in Slit 1 
(upper right). The b/g quasar spectra are located in both Slit 2 (lower left) and Slit 3 (lower right) at the 
same spatial offset of 176 kpc. The spectroscopic observations demonstrate the extreme kinematics of 
the system: AGN 1 has a velocity of +1300 ± 400 km s _1 relative to the f/g quasar and the Lycc emission 
in the nebula exhibits motions ranging from —800km s” 1 (at ~ 100 kpc offset in Slit 3, lower right) to 
+2500km 1 (at ~ 100 kpc in Slit 1, upper right). A 3 X 3 pixel boxcar smoothing, which corresponds 

120 km s -1 X 0.8", has been applied to the spectra. In each two-dimensional spectrum, the zero velocity 
corresponds to the systemic redshift of the f/g quasar. The color bars indicate the flux levels in units of 
erg s -1 cm -2 arcsec -2 A -1 . 9 
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Fig. 4: Absorption line spectrum of cool gas in SDSSJ0841+3921 . Spectrum of the absorbing 
gas detected in the b/g quasar sightline at an impact parameter of 176 kpc from the f/g quasar. The 
gas shows strong Hi and low-ionization state metal absoiption, offset by ~ 650km s -1 from the f/g 
quasar’s systemic redshift. The CII absoiption in particular exhibits a significant tail to velocities as 
large as ~ 1000 km s -1 , providing evidence for extreme gas kinematics. We modeled the strong HI 
Lya absorption with a Voigt profile (blue curve with grey band indicating uncertainty) and estimate a 
column density log Api = 19.2 ± 0.3. The strong low and intermediate ion absorption (Sill, CII, Silll) 
and correspondingly weak high-ion absorption (CIV, SilV) indicate that the gas is not highly ionized, 
and our photoionization modeling (45) implies log 10 rrni = —1.2 ± 0.3 or log 10 Ah = 20.4 ± 0.4. We 
estimate a conservative lower-limit on the gas metallicity to be 1/10 of the solar value. 
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1 Optical Observations 

1.1 Discovery Spectra 

The source SDSSJ084158.47+392121.0 (f/g quasar) was targeted as a quasar by the Sloan Digi¬ 
tal Sky Survey (SDSS) through their color-selection algorithms and observed spectroscopically 
in the standard survey (52). It has a cataloged emission redshift of z — 2.046. Our analysis of 
the SDSS imaging data revealed a second, neighboring source at SDSSJ084159.26+392140.0 
(b/g quasar) with colors characteristic of a z ~ 2 quasar, characterizing this system as a can¬ 
didate quasar pair. In the course of our ongoing spectroscopic campaign to discover close 
quasar pairs at z > 2 (6, 53), we confirmed this source to have z = 2.214 implying a pro¬ 
jected quasar pair with a physical separation of R± = 176 kpc at the f/g quasar redshift. 
SDSSJ084159.26+392140.0 was also observed by the SDSS-III survey in the Baryonic Oscil¬ 
lating Spectroscopic Survey campaign at a spectral resolution of R & 2000 and with wavelength 
coverage A ~ 3600 — 10, 000A (54). 

On UT 2007 Jan 18, we observed the quasar pair with the Low Resolution Imaging Spec¬ 
trograph (LRIS) (55). These data were taken to study the intergalactic medium probed by the 
quasar pair, to examine the HI gas associated with the circumgalactic medium of the f/g quasar, 
and to search for fluorescent emission associated with the f/g quasar. The latter analyses of 
these data have been presented in previous works (22, 26-28). Summarizing the observations, 
we used LRIS in multi-slit mode with a custom designed slitmask which allowed placement 
of one slit on the known quasars and other slits on additional quasar candidates in the field. 
Specifically, a slit was placed at the position angle PA=25.8° between the f/g and b/g quasars, 
allowing them to be observed simultaneously. LRIS is a double spectrograph with two arms 
giving simultaneous coverage of the near-UV and optical. We used the D460 dichroic with the 
1200 lines mm -1 grism blazed at 3400 A on the blue side, resulting in wavelength coverage 
of « 3250 — 4300 A. The dispersion of this grism is 0.50 A per pixel and our 1" slits give a 
resolution with full width at half maximum (FWHM) FWHM~ 160km s~ 1 . On the red side we 
used the R1200/5000 (covering ~ 4700 —6000A) and R300/5000 (covering ~ 4700 — 10, 000A) 
gratings having a FWHM of « 100km s -1 and 400km s’ 1 respectively. 
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The science frames were complemented by a series of calibration images: arc lamp, dome 
flat, twilight sky, and standard star spectra with the same instrument configurations. All of these 
exposures were reduced using the LowRedux (http://www.ucolick.org/~xavier/LowRedux/) 
pipeline which bias subtracts and flat fields the images, corrects for non-uniform illumination, 
derives a wavelength solution, performs sky subtraction, optimally extracts the sources, and 
fluxes the resultant spectra. The ID spectra are corrected for instrument flexure and shifted to a 
heliocentric, vacuum-corrected system. 

Using custom software, we also coadded the 2D spectral images to search for diffuse Lyct 
emission surrounding the f/g quasar (22). This software enables us to model and subtract the 
spectral PSF of sources in the 2D spectra. Extended Lycc emission will be manifest as residual 
flux in our 2D sky-and-PSF-subtracted images which is inconsistent with being noise. To vi¬ 
sually assess the statistical significance of any putative emission feature, we define a x image 
Xsky+PSF = (DATA — SKY — OBJECTS) /a. If our model is an accurate description of the 
data, the distribution of pixel values in the y s k y +PSF should be a Gaussian with unit variance. 
The middle row of the images in Fig. SI shows this quantity for SDSSJ 0841+3921 for each 
slit orientation. The lower row of these images shows y s k y = (DATA — SKY) /a. The upper 
row shows a smoothed version of the middle row, helpful for identifying extended emission. 
Specifically, the smoothed images are given by yy mth = CONx - objects] , w ^ ere 

the CONVOL operation denotes smoothing of the stacked images with a symmetric Gaus¬ 
sian kernel (same spatial and spectral widths) with FWHM=235 km s -1 (dispersion cr smt h = 
100 km s -1 ), which is 5.7 pixels, or 1.4 times the spectral resolution element, and corresponds 
to FWHM = 1.5" spatially. The operation CONVOL 2 represents convolution with the square 
of the smoothing kernel. With this definition of yy mtll the distribution of the pixel values in the 
smoothed image will still obey Gaussian statistics, although they are of course correlated, and 
hence not independent. 

Sky and object PSF-subtracted y-maps for all of the slit-orientations that we used to char¬ 
acterize extended emission in SDSSJ 0841+3921 are shown in Fig. SI. Fig. 3 of the Main text, 
shows just the smoothed sky-subtracted images with a color-map chosen to accentuate the faint 
extended emission and a slightly different smoothing, namely a 3 x 3 pixel boxcar smoothing, 
which corresponds 120 km s 1 x 0. 8". The y-maps in Fig. SI enable the reader to objectively 
assess the statistical significance of all emission features in the unsmoothed data. Note that in 
all of these maps, only the PSF model of the f/g and b/g quasars have been subtracted, whereas 
neither AGN1 nor AGN2 were removed from the other slit-orientations (see Fig. 1 in main text). 
Following the calibration procedure described in (22), we deduce the following spectroscopic 
surface brightness limits (la) for the Fya line at the f/g quasar redshift z = 2.0412: Slitl, 
SBi ct = 2.2 x 10 -18 ergs -1 cm -2 arcsec -2 ; Slit2, SBi ct = 2.3 x 10 -18 ergs -1 cm -2 arcsec -2 ; 
and Slit3, SB i a = 4.8 x 10 18 ergs 1 cm 2 arcsec 2 , where these SBs are computed in win¬ 
dows of 700 km s -1 x 1.0", which corresponds to an aperture of 700 km s -1 x 1.0 arcsec 2 on the 
sky because we always used a 1.0" slit. The depth that we attain in our Fya spectroscopy is com¬ 
parable to that achieved by our narrow-band imaging SBi ct = 1.7 x 10 -18 ergs -1 cm -2 arcsec -5 
(see next section). 
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The original observation from the Quasars Probing Quasars (QPQ) (22) survey corresponds 
to slit-orientation Slit2 (see Fig. SI and Fig. 3 in main text). Our initial visual inspection 
of the Lya emission map revealed a bridge of Lya emission along the slit connecting the 
quasar pair. Although the SB varies with position, this structure has a characteristic SB Lya ~ 
10 -1 ' ergs -1 cm -2 arcsec -2 and is detected at ~ 4 — 5a at location R± ~ +100 kpc between 
the two quasars. The emission extends along the slit from —80 kpc below the f/g quasar to 
+150 kpc between the two quasars, and possibly to +250 kpc. With an end-to-end size of 
~ 250 kpc, this represented one of the largest Lyo: emission nebula ever detected, and thus 
motivated the additional imaging and spectroscopic observations, which are the focus of this 
manuscript. 

1.2 Narrow Band Imaging 

We purchased a custom-designed narrow-band (NB) filter from Andover Corporation, sized 
to fit within the grism holder of the Keck/LRISb camera. The filter was tuned to A cen ter = 
3700A, and designed with a narrow band-pass FWIIM A = 33A to minimize sky background 
while maintaining throughput. On UT 12 November 2012, we imaged the ~ 5' x 7' field-of- 
view surrounding SDSSJ0841+3921, offset to place the quasar pair on the CCD with highest 
quantum efficiency. We observed for a total of 3.2 hours in a series of dithered, 1280s exposures. 
Conditions were clear with sub-arcsecond atmospheric seeing. In parallel, we obtained 3hrs of 
broad-band V images with the LRISr camera. The instrument was configured with the D460 
dichroic and the detectors of the blue camera were binned 2x2 to minimize read noise. 

The images were reduced using standard routines within the IRAF reduction software pack¬ 
age. This includes bias subtraction, flat fielding and an illumination correction. A combination 
of twilight sky flats and unregistered science frames were used to produce flat-field images and 
illumination corrections in each band. The individual frames were sequentially registered to the 
SDSS-DR7 catalog using the SExtractor ( 56) and SCAMP (57) packages. The RMS uncertainty 
in the astrometry of our registered images is approximately 0.2". Finally, the corrected frames 
from each band were average-combined using the SWarp package (55). 

We calibrated the photometry of our images as follows. At the beginning and end of the 
night, we observed two spectrophotometric stars (G191b2b and Feige34) with the NB filter 
under clear conditions. Neither star has significant spectral features in the relevant wavelength 
range covered by our NB filter. For the broad-band images, we observed the standard star 
field PG0231+051. To compute the zero-point for the narrow-band images, we compared the 
measured count rates of Feige34 and G191b2b with the expected fluxes estimated by convolving 
the standard star spectra (resolution of lA) (59) with the normalized filter transmission curve. 
We calculate an average, zero-point magnitude of 23.58 mag, with a few percent difference 
between the two stars. With this calibration, we deduce that the la limiting surface brightness 
of our combined NB images is SBi ct = 1.7 x 10 -18 ergs -1 cm -2 arcsec -2 for a 1.0 arcsec 2 
aperture. 

For the broad-band images, we compared the number of counts per second of the five stars in 
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the PG0231+051 -field with their tabulated R-band magnitudes (60). The derived zero-point for 
the five stars are consistent to within a few percent and we adopt the average value: V Z p = 28.07. 
Because the Feige34 and the PG0231+051 fields were observed with an airmass (AMps 1.2) 
similar to our science field, we did not correct the individual images before combination. We 
estimate that the correction would be on the order of few percent. The la limiting point source 
magnitude for our combined broad band images is V = XX. 

To isolate the emission in the Lya line we estimated and then subtracted the continuum 
emission underlying the NB3700 filter. We estimate the continuum using the 17-band (not 
affected by the Lya line) and assuming a flat continuum slope (in frequency), i.e. (3\ = —2. 
This assumption is dictated by our dataset, i.e. we can rely only on a deep 17-band image to 
estimate the continuum, and it follows the prescriptions of previous work (61) which showed 
that /3\ ~ —1.9 for the most luminous sources and higher values for the faintest. Thus, assuming 
a flat continuum slope gives us a conservative estimate of the Lya emission and its equivalent 
width. As the V --band image has slightly worse seeing (~ 1.3"), we convolved the NB image to 
the same seeing before subtracting the continuum. The continuum subtraction has been applied 
using the following formula 


Lya 


NB3700 - 


FWHM N B3700 TrNB3700 y 
FWHMy Try ’ 


( 1 ) 


where Ly a is the final subtracted image, NB3700 is the narrow-band image, 17 is the 17-band 
image, and Tr NB 37 oo an d Try are the peak transmission values for the NB3700 and 17-band 
filters, respectively. The result of this procedure is a Lya only narrow-band image, which is 
shown in the left panel of Fig. 1. 


1.3 Follow-up Spectroscopy 

After the discovery of an extended nebula in our NB imaging, we opted to obtain additional 
spectroscopy of both the diffuse nebula and several compact sources near the projected quasar 
pair during the same observing run. On UT 13 November 2012 in clear observing conditions 
and seeing varying from FWHM = 0.6 — 1.0", we used the LRIS spectrometer with a 1.0" 
longslit, configured with the D460 dichroic, B1200/3400 grism, and R600/7500 grating. We 
obtained two additional orientations of the longslit as illustrated in Fig. 3: (Slit 1) along the line 
connecting source f/g quasar and AGN2 with a PA=227°; (Slit 3) along the line connecting the 
b/g quasar and AGN1 with a PA=343°. We exposed for a total of 4800s and 2400s for Slit 1 
and Slit 3, respectively. The PA of Slit 1 was chosen to be aligned with the f/g quasar and 
a nearby radio source detected in the FIRST survey (see Supp. 2.2), which we learned about 
just prior to the November 2012 Keck observations. This turned out to be an AGN at the same 
redshift (AGN2). The slit orientation also conveniently allowed us to simultaneously observe a 
bright spot in the Lya image, which is also the companion AGN 1. The orientation of Slit 3 was 
chosen to be aligned with the b/g quasar as a position reference and two other bright spots in 
the nebula, one of which was also covered by Slit 1, i.e. AGN1. Coordinates, photometry, and 
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other information about these targets are listed in Table S1. Corresponding calibration frames 
were obtained and these data were reduced with the same procedures described above. The slit 
orientation of our original January 2007 spectroscopy is also indicated in Fig. 3 as Slit 2, which 
simultaneously observed the f/g and b/g quasars at PA=25.8° for an exposure time of 1800s (but 
in better conditions). 

Our Lya image (left panel of Fig. 1) reveals several compact Lyo-cmitting sources (LAEs) 
with rest-frame Wj jya > 20A in the vicinity of the f/g quasar, which are hence very likely to be at 
the same redshift. On UT 2012 Dec 14, we targeted two of these LAEs with the DEIMOS spec¬ 
trometer (62) on the Keck II telescope in partly cloudy conditions. Specifically, we employed 
the 1.0" long-slit mask oriented to cover the object labeled AGN3 in Fig. 1 (see also Table SI), 
and another LAE which we designate as Targetl at RA=08:41:58.8, DEC= +39:21:57.4, 
which lies off of the image in Fig. 3. The instrument was configured with the 600ZD grat¬ 
ing tilted to a central wavelength of 7200A which provides coverage from « 4600 — 9800A 
with a spectral resolution FWFIM^ 235km s -1 . We took two exposures totaling 4500s. These 
data were reduced and extracted with the SPEC2D pipeline (63, 64) and fluxed using a spec¬ 
trum of the spectrophotometric standard star Feige 34 taken that (non-photometric) night with 
the same instrument configuration. Note that given the limited blue sensitivity of DEIMOS, 
these spectra did not cover the Lya emission line at ~ 3700A. Owing to its faint continuum 
magnitude (V = 25.2 mag) the spectrum of Targetl was inconclusive, although its large Lya 
equivalent width I+r jyQ = 97A suggests it is a real LAE. A setup covering Lya is likely neces¬ 
sary to spectroscopically confirm this source. The object AGN3 is an AGN at the same redshift 
as the f/g quasar, as described in the next section. 


2 Discovery of Four AGN 

Our follow-up spectroscopy reveals three additional AGN within 6 < 18" or R± < 150 kpc 
of the f/g quasar, with very similar redshifts. These objects are labeled as AGN1-3 in Fig. 3, 
and their optical spectra are shown in Fig. S2. Relevant information for all five of the AGN 
associated with the nebula, namely the f/g quasar, the three AGN discovered at the same redshift, 
and the b/g quasar are provided in Table S1. 

Motivated by the discovery of these AGN, we force-photometered the SDSS and WISE 
survey images at the coordinates of each AGN, measured from our deep Keck images. This 
was performed using a custom algorithm (65), which simultaneously models the input sources 
at their respective locations, as well as all other nearby sources detected in the SDSS survey 
imaging. This modeling fully accounts for the potentially overlapping PSFs of the sources 
(a significant issue for WISE given its FWHM^ 6" PSF), and thus effectively de-blends the 
photometry. It also utili z es data from multiple visits, thus effectively co-adding all epochs and 
improving the effective depth over that in the published catalogs. Table S2 shows photometry 
for the five AGN. There we list the SDSS ugriz and WISE 1+1,1+2,1+3,1+4 measurements 
determined from our forced photometry procedure, measurements or limits on the Peak 20cm 
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radio flux F 2 o C m from the FIRST survey, the LRIS R-band photometry, the Lya line-flux, and 
the rest-frame Lya equivalent width. 

2.1 AGN1: An Obscured Type-2 Quasar 

The object AGN1 is embedded in a bright ridge of the nebular Lya emission, as shown in our 
Lya image and the 2D spectrum (Fig. 3). Our LRIS spectra have relatively low continuum 
signal-to-noise S/N ~ 1 — 2 ratio, consistent with the faint continuum magnitude of AGN1: 
V = 23.95 ± 0.07. Nevertheless, we detect strong emission in several lines: HellIA1640, 
C ill]A1908, CII]A2328, and a tentative detection of CIVA1549, along with bright Lycc, which 
our narrow-band imaging indicated should be strong. The FWHM, line flux, and rest-frame 
equivalent widths (EWs) of these lines are listed in Table S3. The presence of emission in 
several high-ionization UV lines (i.e. He II and C IV) characteristic of AGN, clearly indicates 
that this source is powered by a hard-ionizing spectrum, rather than star-formation. Although 
high-ionization lines l ik e C IV, He II, and C III] can be formed in starbursts, stellar winds, the 
photospheres of massive stars, and the interstellar medium, the spectra of star-forming galaxies 
typically exhibit much smaller rest-frame EWs (< 2A) (66) in these lines than observed in our 
spectra. Several independent lines of argument suggest that this object is actually a luminous but 
obscured quasar, also referred to as a Type-2 quasar, as we elaborate on below. From the weak 
He nA1640 line, we measure z = 2.05476 which is ~ 300km s -1 lower than the estimates from 
Lya, and C ill] A1908. We estimate an uncertainty in this redshift of 400km s -1 , which arises 
both because of line-centroiding error, and possible systematic uncertainty about the degree to 
which a high-ionization line like HellA1640 traces the systemic frame for a Type-2 AGN. 

According to unified models of AGN (67, 68), orientation is the primary determining fac¬ 
tor governing the ultraviolet/optical appearance of an AGN. In this context, vantage points 
which are not extincted by the obscuring torus observe a spatially unresolved power-law ul¬ 
traviolet and optical continuum believed to emerge from the SMBH accretion disk, and broad 
high-ionization emission lines (FWHM~ 5000 — 20,000km s -1 ) from photoionized gas in 
the so called broad line region at distances ~ 1 pc from the SMBH. However, all vantage 
points, observe narrow emission lines from the spatially extended (~ lOOpc — 10 kpc) nar¬ 
row (FWHM~ 500 — 1500km s -1 ) line region. When the line-of-sight to the central engine is 
blocked by a dusty obscuring torus, only the narrow line region is seen. AGN lacking broad 
emission lines, but which nevertheless exhibit narrow high-ionization lines are classified as 
Type-2 systems, while those with broad-line emission are classified as Type-1. Alternatives to 
the AGN unification viewing angle picture argue instead that Type-Is and Type-2s represent 
different phases of quasar evolution (69), with all quasars passing through an obscured phase 
before outflows expel the obscuring material. 

Characteristics of Type-2 quasars at other wavelengths are a hardened X-ray spectrum, re¬ 
sulting from photoelectric absorption of soft X-rays presumably from the same obscuring torus 
extincting the broad-line region, and strong mid-IR (A ~ 10 /mi) emission, which represents 
ultraviolet/optical energy emitted from the disk, but absorbed, reprocessed, and re-emitted 
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isotropically by the torus. In a subset of radio-loud systems, synchrotron radiation, believed 
to be emitted isotropically from scales much larger than the torus, is also observed. 

Luminous high-redshift radio-loud Type-2 quasars have been studied for decades (see the 
review by (70)) but their radio-quiet counterparts have been much harder to identify. Over the 
past decade, significant progress has been made in identifying sizable samples of low-redshift 
(z < 1) radio-quiet Type-2s from mid-IR (71, 72), X-ray (73), and optical narrow-line emission 
(74). However, to date the vast majority of these Type-2s are at z < 1, and there are still 
only handfuls of bonafied Type-2 quasars known at z ~ 2 (75) with bolometric luminosities 
comparable to the typical SDSS/BOSS Type-1 quasar, L bol > 10 46 ergs -1 

The three lines of argument that we use to establish that AGN1 is a Type-2 quasar are 1) the 
narrow velocity width of its emission lines 2) the line ratios in its spectrum, and 3) its extremely 
red optical-to-mid-IR colors. 

Unlike the three other AGN physically associated with the nebula (f/g quasar, AGN2, and 
AGN3), which exhibit broad emission lines (FWHM~ 5000 — 10,000km s _1 ), Gaussian fits to 
the lines in AGN 1 reveal relatively narrow lines < 1500km s -1 , which are listed in Table S3. At 
low-redshifts (z < 0.3), AGN exhibit a bimodal distribution of Ha emission line widths, which 
corresponds to the Type-1 and Type2 populations. Based on this distribution, the canonical 
dividing line between Type-l/Type-2 quasar classification is FWHM Ha < 1200 km s -1 (76). 
However, it is unclear how to extrapolate this classification to the rest-frame UV lines observed 
at z ~ 2. First, higher redshift sources are intrinsically much more luminous, have more massive 
SMBHs, and hence SMBH-galaxy scaling relations imply that they should be hosted by much 
more massive galaxies. Indeed, for the handful of Type-1 quasars at z ~ 2 for which the 
narrow-line region [O ill] emission line has been observed in the near-IR, larger line-widths 
FWHM= 1300 — 2100km s -1 are indeed observed (77), contrasting with narrower line widths 
observed at lower redshift (76). Furthermore, even in the narrow-line region, one generally 
expects larger widths for Lya and higher ionization potential lines, because photoionization 
modeling generally predicts that emission in these lines is produced at smaller distances from 
the nucleus, where velocities are likely to be larger. Indeed, exactly such trends were observed in 
the spectrum of a Type-2 quasar at z — 3.288 discovered by (78), who measured FWHM Lya = 
1520 ± 30km s _1 , FWHM He ii = 940 ± 140km s -1 , and FWHM cm ] = 1090 ± 140km s -1 , 
which are comparable to our measured values for AGN1. A near-IR spectrum of the Stem et al. 
source reveals narrower widths for the rest-frame optical lines that are typically used to classify 
low-redshift Type-2s: FWHM^ = 170 ± 130km s -1 and FWHMpm] = 430 ± 30km s -1 . 
Note that the unusually hard X-ray spectrum of the Stem et al. source (78) and the implied 
high photoelectric absorbing column Ajf = (4.8 ± 2.1) x 10 23 cm 2 leave little doubt that it is 
a bonafied Type-2 quasar, notwithstanding the fact that Lya and other high-ionization UV lines 
have line widths in excess of FWHM = 1200 km s -1 . We suspect that the canonical rest-frame 
optical line-width classification based on z < 0.3 Type-2 quasars (76) needs to be revised for 
rest-frame UV lines in more luminous Type-2s at z ~ 2. Based on the narrow lines of AGN1 
and their similarity to the Stern et al. source (78), we conclude that it is very likely to be a 
Type-2 quasar. 
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Table S3 presents measurements of the line fluxes and measurements or limits on rest-frame 
equivalent width for the characteristic AGN lines covered by our spectrum. Based on these 
line flux measurements, Fig. S3 shows the line-ratios of AGN 1 in the CIV /He II vs CIv/C III] 
plane (left) as well as the C m]/He II vs C III]/C II] plane (right). These are the standard line 
ratios conventionally discussed in studies that use photoionization modeling to diagnose the 
physical conditions in the narrow-line region (79), although we note that the ratios C III] /He II 
vs Cm]/Cll] are not typically plotted against each other. In Fig. S3, we compare the line 
ratios of AGN1 to other Type-2 AGN compiled from the literature. Specifically, the circles are 
individual measurements of HzRGs (black) from the compilation of (80), and narrow-line X-ray 
sources (cyan) and Seyfert-2s (blue) from the compilation of (81). In addition, triangles indicate 
measurements from the composite spectra of HzRGs (orange) from (70), the composite Type-2 
AGN spectra (purple) from (82), who split their population into two samples above and below 
Lya EW of 63A, and a composite spectrum of mid-IR selected Type-2 AGN (green) from (S3). 
Finally, the stars represent measurements of these line ratios from composite spectra of Type-1 
quasars, based on the analysis of (magenta) (84), (red-magenta) (35), and (blue-magenta) (86). 

We caution that determination of robust emission line fluxes for Type-1 quasars is extremely 
challenging. Generally the broad emission lines, non-trivial emission line shapes, and the fact 
that many of the lines of interest are severely blended, makes the separation of the spectrum 
into line and continuum rather ambiguous. To minimize the impact of noise and variations 
among quasars, one typically analyzes extremely high signal-to-noise ratio composite spectra, 
which also average down quasar-to-quasar variation (84-86). But nevertheless, the resulting 
line fluxes and hence line flux ratios are dependent on the method adopted (see (35)) section 
4.1 for a detailed discussion). These issues are particularly acute for the lines that we consider: 
HellA1640 is heavily blended with the red wing of Civ A1549 and Olll]A1663, AlllA1670, 
and a broad Fell complex, Cm]A1909 is severely blended with AlllA 1857 and SiHlA 1892, 
and Cll]A2326 is blended with a broad Fell complex. These ambiguities are reflected in the 
large variation in the line ratios determined from Type-1 quasar composites in Fig. S3 (stars) 
by different studies. Despite these challenges, Fig. S3 clearly illustrates that the line-ratios for 
AGN 1 are consistent with the Type-2/Seyfert-2 locus, and differ rather significantly from the 
Type-1 line ratios, notwithstanding the somewhat divergent measurements for the latter. We 
thus conclude based on the emission line ratios of AGN 1 that it is very likely to be a Type-2 
AGN. 

Many studies utilizing WISE data have shown that at bright magnitudes W2 ~ 15 mag, 
the WISE Wl — W2 color efficiently selects AGN at z ~ 1 — 3, including obscured objects 
(87-89). This results from the rising roughly power law shape of AGN SEDs in the mid- 
IR (90,91), and the fact that stellar contamination is low, because the WISE W 1 and W2 bands 
cover the Rayleigh-Jeans tail of emission from Galactic stars, resulting in much bluer colors 
which cleanly separate from redder AGN. At fainter magnitudes W2 ~ 17, contamination 
increases because of optically faint, high-redshift elliptical and Sbc galaxies, which tend to also 
have very red Wl — W2 colors. In addition, obscured AGN are also expected to have very 
red optical to mid-IR colors, if the mid-IR traces reprocessed dust emission from a heavily 
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extincted accretion disk emitting in the optical/UV. Indeed, it has been shown that the color- 
cut r — W 2 > 6 effectively separates obscured AGN from the their much bluer unobscured 
counterparts (89, 92, 93). The WISE photometry in Table S2 indicate that AGN1 is detected in 
both W 1 and W2 with W 1 — 16.89 ± 0.08, and W2 = 15.86 ± 0.11, implying a Wl — W2 = 
1.03 ± 0.14, thus consistent with the expected red color Wl — W2 > 0.8 of AGN in the WISE 
bands. If we adopt our LRIS V-band magnitude as a proxy for r (color corrections are minimal), 
then the faint V = 23.95 continuum of AGN1 implies r — W2 = 8.09 ± 0.13, which is also 
consistent with the r — W 2 > 6 criterion for obscured AGN. We thus conclude based on mid-IR 
and mid-IR to optical colors, that AGN1 is consistent with being an obscured AGN. 

2.2 AGN2: A Radio Loud Broad Line AGN 

The second quasar (AGN2) is rather faint V = 23.12 (Table SI), but its LRIS spectrum never¬ 
theless reveals broad Lya, C IVA1549, and C ill] A1909 emission lines, unequivocally indicating 
that it is a broad-line Type-1 AGN (see Fig. S2). Centroiding the broad and relatively low S/N 
C ill] A1908 line, we estimate ^agn 2 = 2.058 with a 700km s -1 uncertainty. Redshift errors for 
rest-frame ultraviolet emission lines are dominated by intrinsic shifts between the emission line 
and the true systemic frame, and we adopt the procedure described in (23) to centroid the broad 
lines, and use the procedure of (94) to determine the redshift errors. 

Motivated, in part, by the previous association of bright Lya nebulae to HzRGs and/or bright 
radio sources, we searched the literature and public catalogs for sources in the field surrounding 
SDSSJ0841+3921. The Faint Images of the Radio Sky at Twenty cm survey (FIRST) (95) 
used the Very Large Array (VLA) to produce a map of the 20 cm (1.4 GHz) sky with a beam 
size of 5. "4 and an RMS sensitivity of about 0.15 mJy beam -1 . The survey covers the same 
10,000 deg 2 sky region covered by the SDSS imaging, has a typical detection threshold of 
1 mJy, and an astrometric accuracy of 0.05". The FIRST catalog reports only a single source 
within 60" of SDSSJ0841+3921 at RA=08:41:58.6 DEC=+39:21:14.7, approximately 6" South 
of the f/g quasar and coincident with AGN2. The FIRST peak and integrated fluxes at 20 cm are 
Fpeak = 12.87 mJy/beam and F int =14.81 mJy respectively, giving a morphological parameter 
0 = log 10 (F int /Fp eak ) = 0.06. This indicates a compact, unresolved source at 1.4 GHz (96). 
Therefore, there is no evidence for a bright, extended radio source (e.g. a radio jet) associated 
with the protocluster system and giant Lya nebula. 

Follow-up VLA observations of this source were performed by (97), having erroneously as¬ 
sociated it with the f/g quasar. The spectral slope measured from 8.4 GHz and 4.9 GHz observa¬ 
tions is relatively steep, af g = —1.21, which is generally interpreted as evidence that the source 
has a large viewing angle (i.e. edge-on) (98). This may indicate that FIRSTJ084158.6+392114.7 
(AGN2) is oriented away from us, although we re-emphasize that it exhibits a Type-1 spectrum. 
These much deeper radio images of the source show no evidence for extended structure or jets, 
providing further evidence that it is not a HzRG (97). 
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2.3 AGN3: A Faint Broad Line AGN 


Similar to AGN2, the source AGN3 is also faint V = 23.09 but exhibits broad emission-lines of 
CIVA 1549, C III] A 1908, and Mg ilA 2798 (see Fig. S2). Unfortunately Lyct was not covered by 
our DEIMOS spectrum, but the source has a large rest-frame equivalent width (VF Lya = 35A) 
indicative of strong emission. Based on the broad lines, we also classify this source as a broad¬ 
line (Type-1) AGN. We determine a redshift of zagn 3 = 2.050 from the C ill] A1908 line with a 
700km s' 1 uncertainty. 


3 Probability of Finding a Quadruple Quasar 

The probability dP of detecting three quasars around a single known quasar (here the f/g quasar) 
can be written as 

dP = n% so dV 2 dV 3 dV 4 [1 

T£(?T 2 ) + --- (6 permutations) 

+ C(' r i 2 ,r’ 23 ,r' 3 i) H- (4 permutations) (2) 

+ (((r* 12 )^ 34 ) + ••• (3 permutations) 

+ v( r lU 2U3U4)] 

where tiqso is the number density of quasars, dV, is the infinitesimal volume element centered 
on the location r* of the 2 th quasar, is the distance between the /th and jth quasar, and £, 
and r) are the two-point, three-point, and four-point correlation functions, respectively (99). 

The total probability P of finding a quadruple quasar system within a maximum radius r max 
is the integral f dP over the three volume elements dV r , for all possible configurations of the f t . 
Our goal in what follows is to obtain an order of magnitude estimate for this total probability. 

First note that on the small scales of interest to us (here r < 200 kpc), many of the terms in 
eqn. ([2]) can be neglected. For the higher order correlation functions ( and ?/, a scaling of the 
form ( oc £ 2 and ?/ oc £ 3 is often assumed. This arises from the fact that for Gaussian initial 
conditions and a scale-free initial power spectrum, this scaling is obeyed to second order in 
Eulerian perturbation theory (100). However, this scaling is not valid for the highly non-linear 
small-scales of interest to us here. On these small scales, we follow the halo model approach 
(101), and write the higher order functions as a sum of terms representing contributions from 
the multiple possible halos. For example the four-point function can be written as the sum 

V = + n 2h + rf h + 27 4h , (3) 

resulting from contributions from one to four halos. As we are concerned with scales com¬ 
parable to the virial radius of the dark matter halo hosting, i.e. the f/g quasar r < 200 kpc, 
the one halo term, which quantifies the four-point correlations when all four points lie in the 
same dark matter halo, will dominate (101). It is easier to calculate this one-halo term of the 
four-point function ?/ lh in Fourier space, where one instead works with the one halo term of the 
trispectrum T lh , and it can be shown that T lh scales as the Fourier transform of the dark matter 
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halo density profile to the fourth power T lh oc p^ M . This scaling holds likewise for r/. Thus for 
quasar 4-point correlations on small-scales, the halo model indicates that we expect 77 oc Pq SO , 
where Pqso is the density profile of quasars in the dark matter halos which host them. 

By a similar line of argument, the small-scale two- and three-point functions scale as £ oc 
Pq SO and (' oc Pq SO , respectively. Or alternatively, (' oc £ 3 / 2 and r/ oc £ 2 . On the proper scales 
of interest to us here r < 200 kpc (r com = 600kpc, comoving), £ 1 and thus the terms 

dominating eqn. ([2]) are the three permutations of terms like £(r 12 )£(r 34 ), and //, which are all 
of comparable order £ 2 . 

Thus in order to compute the total probability of finding three quasars around a known 
quasar, we must integrate the four dominant terms in eqn. Q, over all possible configurations 
of the Tj, which can be realized within a spherical volume, whose radius is set by the largest 
separation in the quadruple quasar r max . These integrals all have a common form and order- 
of-magnitude. For example the first three terms involving pairs of correlation functions can be 
written 


n 


QSO 


£(r)47rr 2 dr 


C(r M )dV :i dV 4 ~ n\ 


^QSO 


£(r)47rr 2 c/r ) V, (4) 


where V = 47r/3rf nax . Thus the final expression for an order-of-magnitude estimate of the 
probability is 

( pr max \ 2 

J 1 £(r )4:7Tr 2 drJ V , (5) 

where the factor of four accounts for the four similar dominant terms of order ~ £ 2 in eqn. (|2j). 

To calculate n qso, we use the quasar luminosity function of (102) assuming an apparent 
magnitude limit V < 24. This limit is highly conservative, as even AGN2 and AGN3 have 
V = 23.1 mag. Although the Type-2 AGN1 has V = 24, its nucleus is obscured from our 
perspective, and this magnitude represents host-galaxy emission or a small-fraction ~ 5% of 
scattered nuclear emission. The corresponding nuclear V -band emission of this AGN is likely 
to be much brighter than V = 24 given that its Lyo: flux, coming from the narrow-line region, is 
so strong, i.e. larger than that of AGN2. To account for the fact that our luminosity function only 
includes unobscured AGN, we divide it by a factor (1 — f Qbs ), where / obs is the obscured fraction 
of AGN. We adopt an obscured fraction of f Q bs = 0.5 consistent with recent determinations 
(103,104). Altogether, we estimate that the comoving number density of faint AGN is 71qso = 
3.8 x 10 -5 cMpc~ 3 . 

We assume that all four AGN reside in a physical structure with physical size r max . This is 
a good assumption, given that: the projected separations of the AGN are small, they reside in a 
larger-scale overdensity of Lycc emitters, the f/g quasar, AGN1, and AGN2 are all embedded in 
the giant Lycc nebula, and the morphology of the nebular emission around the AGN suggests a 
physical association. Moreover, although we argue here that the probability of finding a quadru¬ 
ple AGN is very small, it would be much smaller if the AGN had larger line-of-sight separations, 
in which case one cannot invoke as large of an enhancement due to clustering. Taking the f/g 
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quasar location as the origin, the largest angular separation measured in the quadruple quasar is 
that of AGN3 with 6 = 18", corresponding to to a transverse separation r± = 150 kpc. We thus 
choose r max = 250 kpc, such that the the spherical comoving volume we consider corresponds 
to V = 1.9 cMpc 3 . The probability of finding the three AGN around the f/g quasar at random, 
i.e. in the absence of clustering, is extremely small P ran = (nQso ^) 3 = 3.4 x 10 -13 . 

The small-scale two-point correlation function of quasars was first computed by ( 6 ), and 
later extended to even smaller scales ~ 5 kpc by the gravitational lensing search of (30), who 
found that a power law form £ = (r/r 0 ) -7 with 7 = —2 and r 0 = 5.4 ± 0.3 /r _1 cMpc provides 
a good fit to the data over a large range of scales. Plugging these numbers into eqn. 0, we 
finally arrive at a probability of P = 6.9 x 10 -8 , justifying our order of magnitude estimate of 
P ~ lx 10~ 7 . 

4 Giant Nebulae-LAE Clustering Analysis 

4.1 Constructing the LAE Catalog 

We created an LAE catalog from our images of SDSSJ0841+3921 using SExtractor (105) in 
dual-mode, using the NB image as the reference image. In order to minimize spurious detection 
we varied the parameters DETECT JYIINAREA (minimum detection area) and 
DETECT.TIIRESII (relative detection threshold) creating a large number of LAE candidate 
catalogs. We verify the number of spurious detections for each catalog running SExtractor on 
the “negative” NB image obtained by multiplying the NB image by —1. The ratio between the 
number of sources detected in the NB image and the sum of sources detected in both the NB 
and “negative” image define a “reliability” criterion. We selected DETECT_MINAREA equal 
to 8 pixels and DETECT.THRESH equal to 1.8cr corresponding to a ’’reliability” of 94%. We 
used the same parameters to obtain a catalog of sources for the E-band image. In order to select 
LAE candidates, we estimated the NB — V color from SExtractor isophotal magnitudes and 
applied a color-cut corresponding to a Lya rest-frame equivalent width of W\, yoi > 20A which 
is the standard threshold used in the literature. In particular, we assumed a flat continuum slope 
in frequency to convert the V-band magnitude to a continuum flux at the Lya wavelength. For 
the objects in the catalog without broad-band detection (i.e., with S/N< 3), we determine a 
lower limit on the EW following (61). Finally, we inspected each object by eye and removed 
four spurious sources that were lying in proximity to a bright star. 

The final clean catalog consisted of 61 LAE candidates above a Lya flux of limit of 5.4x 1CT 18 
erg s _1 cm -2 , over our 6 ' x 7.8' imaging field-of-view, and the redshift range z = 2.030 — 2.057 
spanned by the narrow-band filter. We estimated that our catalog is 50% (90%) complete above 
a flux limit of P 50 = 6.7 x 10 -18 erg s _1 cm -2 (P 90 = 7.4 x 10 -18 erg s _1 cm -2 ). This catalog 
of LAEs is presented in Table S4. In the following clustering analysis, we consider only those 
sources with rest-frame Wj Aya > 20A and a luminosity log 10 L^a > 42.1 which at z = 2.04 
corresponds to a flux limit of 4.0 x 10 — 1 ' erg s _1 cm -2 . Given that this flux level is a factor 
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of five higher than the 90% completeness flux limit of our catalog, we are certain that we are 
complete to such sources. The f/g quasar lies at a distance of 2.08' from the edge of our imag¬ 
ing field, and for simplicity, our clustering analysis considers only the 10 sources within this 
radius which corresponds to a comoving impact parameter of 2.2 h~ l cMpc. Of these, 3/10 are 
the spectroscopically confirmed AGN1-3. The sources included in our clustering analysis are 
denoted by an asterisk in Table S4. 

4.2 Clustering Analysis of HzRGs and LABs 

We now describe how we used a compilation of LAEs around HzRGs and LABs to estimate 
the giant nebulae-LAE cross correlation function. The vast majority of z ~ 2 — 3 protoclusters 
in the literature have been identified and studied via overdensities of LAEs over comparable 
fields of view 6' x 7.8' as our LRIS observations. For the HzRGs, we used a catalog of LAE 
positions around HzRGs from the survey of (13). Specifically, we focus on six HzRGs from 
the Venemans et al. study in the redshift range z = 2.06 — 3.13, which are approximately 
co-eval with SDSSJ0841+3921 at z — 2.0412. In addition to the six HzRGs from Venemans 
et al. (13), we also include two LABs at z ~ 2 — 3 from the literature whose environments 
have been surveyed for LAEs. Nestor et al. (106) conducted a narrow band imaging survey 
of the well known SSA22 protocluster field at z — 3.10. Their Keck/LRIS observations were 
centered on LABI, whereas LAB2, the other known LAB in this field, resides at the edge of 
their imaging field. For this reason we only consider the overdensity around LAB 1. Prescott et 
al. (19) conducted an intermediate-band imaging survey for LAEs around an LAB at 2 = 2.66 
(LABd05) discovered by Dey et al. (16), finding a significant overdensity, and argued that the 
LAB resides in a protocluster. The eight objects used in our clustering analysis are listed in 
Table S5. 

Given that the Venemans et al. sample (13) is the largest compilation of LAEs around 
HzRGs, and that, to our knowlege, the two LABs we consider are the only ones whose envi¬ 
ronments have been characterized using LAEs to a depth comparable to our observations, this 
sample of eight HzRGs/LABs is the only dataset currently avaailable for conducting our clus¬ 
tering analysis and we believe that they comprise a fairly representative sample. Venemans et 
al. applied standard HzRG selection criteria to define their sample, and they published results 
for all HzRGs, not just those that hosted dramatic overdensities. Although LABI and LABd05 
are larger and brighter than more typical LABs, they are comparable in size and luminosity to 
the HzRGs and to SDSSJ0841+3921, making for a fair comparison. One cause for concern is 
that the distribution of LAEs around the LABs were published specifically because these two 
systems resided in overdensities, and thus there may be a publication bias against LABs resid¬ 
ing in lower density environments. Given that LABS only comprise a fourth of our clustering 
sample, we do not believe that this significantly biases our results, but repeating our clustering 
analysis for a larger and more uniformly selected sample of LABs would clearly be desirable. 

Note that the HzRG MRC 1138—2262 included in our analysis is the famous Spider Web 
Galaxy protocluster, which has been the subject of extensive multi-wavelength follow-up, and 
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is one of the most dramatic protocluster systems known. Venemans et al. (13) argued that all of 
the HzRGs in Table S5 reside in overdensities of LAEs, with the exception of MRC 2048—272 
and TN J20093040, for which the abundance of LAEs was found to be consistent with the field 
value, albeit with large error bars. We nevertheless include these two sources in our clustering 
analysis for several reasons. Our objective for the clustering analysis is to shed light on the con¬ 
nection between active SMBHs, large-scale Lya nebulae, and protocluster environments. Given 
that it has been argued that the HzRGs as a population trace extremely overdense protocluster 
environments, and given that MRC 2048272 and TN J20093040 are both powerful HzRGs with 
large Lya nebulae, excluding them from the analysis would be arbitrary, and would bias our 
correlation function high. Furthermore, the background number density of LAEs used by Vene¬ 
mans et al. (13) is rather outdated, the overdensity estimates have large error bars for individual 
systems, and the luminosity limits that they adopted for the LAEs in each field considered 
were heterogeneous. These factors complicate the interpretation, particularly if the clustering 
of LAEs is luminosity dependent. Hence our goal is to conduct a homogeneous analysis of 
the clustering of LAEs around the HzRGs/LABs which are above a uniform Lya luminosity 
log 10 Lhya > 42.1, without cherry picking specific objects. We adopt this luminosity limit for 
two reasons. First, the NB imaging observations for all eight objects listed in Table S5 are com¬ 
plete for identifying LAEs above this luminosity, mitigating uncertainties due to incomplete¬ 
ness. Second, this limit corresponds to the faintest luminosity for which the the background 
number density of LAEs is well characterized (32); working fainter would thus require a signif¬ 
icant extrapolation of the luminosity function, which would add systematic uncertainties to our 
results. For all of the objects analyzed (see Table S5), the LAE catalogs available are complete 
above the Lya equivalent width limit of WL ya > 20A (the same value used to define our LAE 
catalog around SDSSJ0841+3921 ), with the exception of the Prescott et al. observations of 
the LAB LABd05. These observations used an intermediate band filter, and thus adopt a rest- 
frame W\, ya > 40A. We consistently account for this difference below when we compute the 
background number density of LAEs. 

Given the small size of this HzRG/LAB dataset, we use an unbinned maximum likeli¬ 
hood estimator to determine the parameters r 0 and 7 , following the procedure outlined in 
(107,108). Specifically, we assume that the cross-correlation function between the giant nebu¬ 
lae (HzRGs/LABs) and LAEs obeys a power law form £(r) = (r/r 0 ) -7 . We use this correlation 
function to calculate the expected number of LAEs within a comoving cylindrical volume with 
transverse separation from R to R + dR, and half-height A Z, which is set by the width of 
the given narrow band filter A Z = FW 11X1/ (2a 11 (z)), where FWHM = cAA/A is the full 
width half maximum of the filter in velocity units, a = 1/(1 + z) is the scale factor, H(z ) is 
the Hubble expansion rate, and division by aH(z ) converts velocities to comoving units. We 
imagine dividing the transverse separation R into a set of infinitesimal bins of width dR, such 
that each bin can contain only one or zero LAEs. Under the assumption that the clustering of 
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LAEs around the giant nebulae is a Poisson process, we can write the likelihood function as 


'Nlae 


*= n 


e 



( 6 ) 




where /i = 47rf?AZn L AE[l + h(R)]dR is the probability of finding a pair in the interval d,R, the 
index i runs over the dR bins containing the N L ae HzRG/LAB-LAE pairs in the sample, and 
the index j runs over all the bins for which there are no pairs. Here tilae(> Lhya) is the number 
density of LAEs brighter than the survey limit L^ ya , and h(Rj is the cross-correlation function 
averaged over the filter width 



(7) 


Taking the natural logarithm of the likelihood, we have 



( 8 ) 


where we have dropped all additive terms which are independent of our model parameters 
(r 0 ,7), and [ /f. tr , in ,7T. m;1X J is the range of comoving scales over which we search for HzRG/LAB- 
LAE pairs. In the following clustering analysis, we adopt R m in = 50 h~ 1 ckpc to minimize 
confusion with LAEs embedded in the respective nebulae, and R max = 2.5 h~ x cMpc, which 
is well matched to the maximum impact parameter probed R = 2.2 h~ l cMpc around the f/g 
quasar. Given that the HzRGs and LABs in Table S5 all have different filter widths, as well as 
different VE Lya limits and redshifts (which changes the background density tilaeX the likelihood 
in eqn. <[8]) is specific to a single source. However, the full likelihood of the entire dataset given 
the model parameters, can be determined by simply summing log-likelihoods (multiplying the 
likelihoods) over all of the HzRGs/LABs, each of which has the form of eqn. ([8]), 

Our clustering analysis relies on a determination of the field LAE luminosity function tilae- 
We use the Ciardullo et al. (32) Schechter function fits to the luminosity function from their own 
LAE survey at z = 3.1, and the Ciardullo et al. fits to the LAE survey data of (109) at z = 2.1. 
Both of these surveys are highly complete down to our limiting luminosity of log 10 L hya > 42.1 
for LAEs with equivalent width W^a > 20A. Ciardullo et al. also estimated the rest-frame 
equivalent width distribution, and found that LAEs follow an exponential distribution with rest- 
frame scale lengths of W 0 = 50A and W 0 = 64A, at z = 2.1 and z — 3.1 respectively. 
Hence, for LABd05 which employed an W limit > 40A, we take the luminosity function of 
LAEs to be n LAE (> W hmit , > L hya ) = exp [-(IE limit - 20A)/IE 0 ]^lae(> 20A, > L LyQ ), 
where tilae(> 20A, > L^ ya ) are the Ciardullo et al. Schechter-function fits. To determine 
the luminosity function at intermediate redshifts between z = 2.1 and z = 3.1, we linearly 
interpolate. 
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In Fig. S4 we show confidence regions in the r 0 — 7 plane for the likelihood in eqn. ([ 8 ]). The 
maximum likelihood value is r 0 = 29.3 h~ l cMpc and 7 = 1.5, however a strong degeneracy 
exists between r 0 and 7 , such that neither is individually well-determined. The low signal-to- 
noise ratio of the clustering measurement is partly to blame for this degeneracy, however it is 
also a consequence of the functional form adopted for the correlation function £(r) = (r/ro) -7 , 
for which the amplitude and slope of the correlation function are not independent parameters. 

Note that errors on the parameters r 0 and 7 given by our maximum likelihood estimator 
in eqn. ([ 8 ]) and illustrated in Fig. S4 underestimate the true errors. This is because our likeli¬ 
hood explicitly assumes that the positions of the LAEs are completely independent draws from 
a Poisson process. This approach ignores correlations between the LAE positions (which are 
generically present for hierarchical clustering) and fluctuations due to cosmic variance. We 
opted for this unbinned estimator, because the traditional method of computing the correlation 
function in impact parameter bins, and fitting a power-law will, for such small samples, yield 
results which depend sensitively on the choice of binning, unless the covariance of the bins is 
included in the fit. The correct approach would be to fit the binned correlation function using the 
correct covariance matrix, whose off-diagonal elements would reflect the correlations between 
LAEs, and whose amplitude would include fluctuations due to sample variance. Given our lim¬ 
ited sample size of eight HzRGs/LABs, the total number of LAEs in our clustering analysis is 
only 75, and there is no hope of computing the covariance from the data itself. In principle one 
could adopt a model for the covariance based on IV-body simulations, but this would require 
a full model for how LAEs populate dark matter halos, and would be sensitive to the assumed 
masses of the protocluster halos. Such a detailed analysis is beyond the scope of the present 
work. As a compromise, we use the unbinned maximum likelihood estimator in eqn. © to 
determine our best fit model parameters, and we obtain estimates of the errors on parameters 
via a bootstrap technique. Specifically, we generate 1000 realizations of mock datasets, where 
a random sampling of the eight HzRGs/LABs in Table S5 are selected with replacement. Given 
that we are measuring a cross-correlation of LAEs around eight distinct objects, only the LAE 
positions arising from the same source are correlated, whereas the others are completely in¬ 
dependent because of the large respective distances between the HzRGs/LABs. As such, for 
each of these 1000 realizations, we also bootstrap resample the positions of the LAEs around 
each object (with replacement), which encapsulates fluctuations due to correlations in the LAE 
positions and Poisson counting fluctuations. Lor each bootstrap realizationn of our dataset, we 
compute the parameters r 0 and 7 via our maximum likelihood procedure, and determine errors 
on these parameters from the resulting distributions. 

Based on this procedure, we define the lex confidence region from the 16th and 84th per¬ 
centiles of the bootstrapped parameter distributions. We find a range 0.36 < 7 < 1.98 for 7 , 
and 15.8h _1 cMpc < r 0 < 9404 h -1 cMpc. The correlation length r 0 is poorly constrained, 
because for flat correlation function slopes 7 < 1 , the clustering projected along the line-of- 
sight is very insensitive to r 0 , which again simply reflects the fact that the amplitude and slope 
of the correlation function are not independent parameters. This is the same behavior seen in 
the contours in Lig. S4. Lixing the slope 7 to its maximum likelihood value of 7 = 1.5, we 
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determine a la confidence region for the cross-correlation length of r 0 = 29.3 ± 4.9 hr 1 cMpc. 

Although we do not fit the binned correlation function to determine model parameters, it is 
instructive to compute a binned correlation function to visualize the data, our model fits, and 
the errors on these fits. As such, we define a dimensionless correlation function as 


X(Rii Ri+1 ) 


A/proto 

E -^Vproto 

ij 


(PG)ij 
(PRhj 


- 1 , 


(9) 


where (PG)ij is the number of real HzRG/LAB-LAE pairs in the impact parameter bin \R t , R i+ 1 ] 
around the j th object. The quantity (PR)ij = U],ae(Rj, > Ibiimitj)^' i- s the expected random 
number of pairs, which is simply the product of the background number density of LAEs, mlae, 
and the cylindrical volume, 14, , defined by the bin [ R,, R i+ x ] and the half-height AZj of the NB 
filter used to survey the jth object. We compute errors on the binned correlation function via 
an analogous bootstrap method, whereby x(R-n Ri+i) is recomputed from 1000 bootstrap re¬ 
samplings of the sample of eight HzRGs/LABs while simultaneously bootstrap resampling the 
LAE separations from each object. 

In Fig. S5, we show the binned correlation function calculated via this procedure, with la 
error bars, determined by taking the 16th and 84th percentiles of the distribution of x(Ri, Ri+i) 
resulting from the bootstrap samples. For any value of the model parameters r 0 and 7 , we can 
compute the predicted value of y R i+ \) by evaluating the integrals 


rAZj rRi+ 1 

(PG)ij = n LAE ( Zj , > Wi imitJ )I4 / / [1 + aVR 2 + Z 2 )}2 -kR dRdZ, (10) 


and then evaluating the sum in eqn. fii. The uncertainty on the predicted clustering level, 
arising from the uncertainty on our model parameters, can thus be determined by evaluating the 
16th to 84th percentile confidence region of the predicted x(Ri, R;+\ ), using the distribution of 
r 0 and 7 from the bootstrap resampling of our model fits. The resulting lcr uncertainty on the 
predicted correlation function is shown as the gray shaded region in Fig. S5, with the maximum 
likelihood value r 0 = 29.3 h _1 cMpc and 7 = 1.5 shown by the red curve. We see that our 
maximum likelihood model fit and bootstrap errors are consistent with the binned correlation 
function estimates and their respective errors. 

Given the results of the clustering analysis described in this sub-section, we generate Fig. 2 
of the main text as follows. Using our catalog of LAEs around the f/g quasar with log 10 L Lya > 
42.1 (Table S4), and the LAE luminosity function prediction for the background number density 
^lae, we can compute the cumulative overdensity profile 5(< R) = 1 + £(-R min , R) = ( PG)/ 
(PR ). The overdensity profile predicted by our maximum likelihood determination of r 0 and 


7 for the giant nebulae-LAE correlation function are determined by evaluating eqn. (10). This 
is shown as the red curve in Fig. 2. Similarly we can determine the lcr error on the predicted 
d(< R) by evaluating it for all of our bootstrap samples of r 0 and 7 and taking the 16th-84th 
percentile region of the resulting d(< R) distribution. This region is shown as the shaded grey 
area in Fig. 2. For comparison we show the expected overdensity (projected over our NB filter) 
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of Lyman break galaxies (LBGs) around radio-quiet quasars, based on the recent measurements 
of quasar-LBG clustering at z ~ 2.7 by Trainor et al. (JO), who found a cross-correlation 
length of r 0 = 7.3 ± 1.3 /r _1 cMpc, assuming a fixed slope of 7 = 1.5, similar to the best-fit 
slope deduced for our giant nebulae-LAE correlation. The blue dotted lines indicate the la 
error quoted by Trainor et al. (10) on the cross-correlation length with the slope held fixed. 
Comparing to the Trainor et al. (10) measurements may seem questionable, given that: a) this 
is the cross-correlation with LBGs, whereas the clustering analysis in Fig. 2 is concerned with 
LAEs; b) the quasars studied by Trainor et al. were hyper-luminous. However, both of these 
differences are likely to produce small effects on the cross-correlation strength because a) the 
clustering of LAEs and LBGs are comparable at z ~ 2.7, and, assuming hierarchical clustering, 
the cross-correlation strength is only sensitive to the square-root of any differences in clustering 
strength; and b) the Trainor et al. (10) results are consistent, within the errors, with the quasar- 
LBG cross-correlation results of (110), who studied much fainter quasars. Indeed, both (110) 
and Trainor et al. (10) argue that the clustering of LBGs around quasars is independent of 
luminosity, assuaging concerns about comparing to hyper-luminous quasar clustering. 


5 Near-IR Observations: The Systemic Redshift of the f/g 
Quasar 

Quasar redshifts determined from the rest-frame ultraviolet emission lines (redshifted into the 
optical at z ~ 2 ) can differ by up to one thousand kilometers per second from the systemic 
frame, because of the complex motions of material in the broad line regions of quasars (94,111). 
To achieve higher precision, one may analyze spectra of the narrow [O in] emission lines emerg¬ 
ing from the narrow-line region, but at z ~ 2, this requires near-IR spectroscopy. On UT 2007 
April 24, we obtained a near-IR spectrum of the f/g quasar (and also the b/g quasar) with the 
Near InfraRed Imager and Spectrometer (NIRI) on the Gemini North telescope. These data 
were reduced and calibrated with the LowRedux package, which performs wavelength calibra¬ 
tion using the atmospheric sky lines recorded in each spectrum. With the chosen configuration, 
the spectra span 1.40 — 1.94//m in the H-band. 

Fig. S 6 presents the NIRI spectra, fluxed with telluric standard stars taken that night. Each 
quasar shows the detection of strong [O ill] emission and modest H (3 emission. The [O ill] A5007 
lines were centroided using a flux-weighted line-centering algorithm with pixel values weighted 
by a Gaussian kernel, and this algorithm was iterated until the line center converged to within 
a specified tolerance (see (25) for more details). A small offset of 27km s -1 was also added 
to the redshifts determined in this manner to account for average blueshift of the [Olll]A5007 
from systemic (111). We measured wavelength centroids of A 0 b s = 1.5231 ± 0.0001/rm and 
1.6095 ± 0.00006// 111 for the f/g and b/g quasars, respectively. This gives systemic redshifts 
% = 2.0412 for the f/g quasar and 2.2138 for the b/g quasar. We estimate an uncertainty in 
this redshift of 50 km s' 1 based on error in our wavelength calibration, line centering, and the 
estimated uncertainty when one uses [O ill] to assess the systemic redshift of a quasar ( 111 ). 
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6 Kinematics 


In the Main text, we described the complex kinematic motions of the quadruple AGN system 
and nebula as measured with slit spectroscopy (optical and near-IR; Figs. 3,S6) and absorption¬ 
line spectroscopy of a b/g quasar (Fig. 4). The AGN exhibit line-of-sight velocity differences 
of at least 8v = 1300km s _1 , requiring extreme gravitational motions. Note that these large 
velocities cannot be explained by Hubble expansion. The miniscule probability of finding a 
quadruple quasar system in the absence of clustering P ~ 10 -13 (see Supp.[|]), and the physical 
association between the AGN and giant nebula, demand that the four AGN reside in a real 
collapsed structure, and thus the relative motions of the AGN must result from gravitational 
motions in a collapsed structure. 

The velocities of the nebular Lya emission are also extreme, exhibiting line-of-sight veloc¬ 
ities relative to the systemic redshift of f/g quasar ranging from 5v = —800 to +2500 km s -1 . 
These also trace the gravitational potential of the protocluster but may be influenced by other 
physical processes (e.g. kinematic or radiative feedback from the AGN). Finally, the absorp¬ 
tion line kinematics of cool metal-enriched gas in the protocluster halo, measured from the b/g 
quasar spectrum, show strong absorption at ~ +650 km s -1 with a significant tail to velocities 
as large as ~ 1000 km s -1 . 

We further emphasize two aspects here. First, although the nebular emission exhibits large, 
relative velocity offsets, the internal motions are relatively quiescent and generally unresolved 
at the Keck/LRIS spectral resolution (FWHM « 160km s -1 ). There is no evidence for the 
“double-peaked” emission characteristic of resonantly trapped Lya. This implies that the gas 
within the nebula has modest opacity to Lya photons and that the observed radiation is not 
predominantly scattered photons. These data also indicate that the gas motions are not strictly 
random, but are rather organized into coherent flows, i.e. the velocity offsets observed greatly 
exceed the measured velocity dispersions. Similar coherent motions have been previously ob¬ 
served in giant Lya nebulae in HzRGs ( 112-115 ) and LABs (116). 

Second, consider the gas kinematics along the single Slit 2 (Figs. 3and SI). This slit covers 
the nebula from near the f/g quasar, where one observes a negligible velocity offset (8v ~ 
0km s -1 ), to the b/g quasar at « 200 kpc separation. As one travels along the slit, the velocity 
offset increases to ~ +500km s -1 as one approaches the b/g quasar, where the absorption-line 
spectrum shows a velocity offset of 8v ~ +650km s” 1 (Fig. 4). Altogether, the gas along 
Slit 2 exhibits a velocity gradient of at least 500km s -1 across 200kpc (projected). These data 
establish the presence of large-scale, high-velocity flows of cool gas within this protocluster 
halo. 
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7 Estimates for the Abundance of Giant Lya Nebulae Around 
Quasars 

The protocluster SDSSJ0841+3921 was discovered because of its large-scale Lya emission 
properties, and was selected from a well-defined survey for extended Lya emission around 29 
quasars (22). We can use the statistics of this survey to estimate the prevalence of giant Lya 
nebulae around quasars. Insofar as there appears to be a physical connection between giant Lya 
nebulae and protoclusters, this also provides insight into the prevalence of protoclusters around 
quasars. We also present an estimate for the frequency of occurrence of giant Lya nebulae from 
an independent narrow-band imaging survey of quasars conducted by our group, and also dis¬ 
cuss the results of similar narrow-band imaging surveys performed by others. Although much 
past work has been dedicated to characterizing Lya nebulae around radio-loud quasar sam¬ 
ples (11 7), it is only recently that sensitive observations have been undertaken to characterize 
the abundance of nebulae around predominantly radio-quiet quasars. 

Our slit-spectroscopic observations allow us to estimate the effective covering factor of ex¬ 
tended Lya emission around the typical quasars (L bo i ~ 10 46 ergs' ') we surveyed. Follow¬ 
ing the discussion in § 4.3 of (22; see Fig. 4 of that paper), large ~ 70 kpc scale nebulae 
can be detected down to 2 x SB !rT , where SB lfJ are the la surface brightness limits quoted 
in Table 1 of (22). In that sample a nebula of comparable brightness to SDSSJ 0841+3921 
SB Lyo , ~ 10 - 17 ergs -1 cm -2 arcsec -2 could have been detected at > 2 a significance in 23/29 
cases (see Table 1 of (22)). Based on this, we can crudely determine the covering factor of 
such filamentary emission as being the ratio of the area of the filament covered by our origi¬ 
nal slit observation of SDSSJ0841+3921 (Slit2 in Figs. 3and SI), to the total area surveyed at 
sensitivity sufficient to detect a SBL ya > 10 -1 ' ergs -1 cm -2 arcsec -2 nebula. 

We consider an annular region around a z = 2 quasar with an inner radius of R ±An \ n = 
50 kpc (5.8") and outer radius of R L , max = 150 kpc (17.3"). The inner radius is chosen to ex¬ 
clude the small-scale fuzz at R± < 50 kpc, which is detected in a large fraction of quasars (22), 
from the larger scale filamentary emission. The outer radius corresponds to the maximum extent 
of the (high significance) emission detected in our discovery spectrum of SDSSJ 0841+3921 
(see Slit2 in Fig. 3panel c, and the middle column of Fig. SI), and is about the virial radius of 
the dark matter halos M vir = 10 12,5 M 0 at z = 2 which host quasars (11). We also observe 
a drop-off in the covering factor of optically thick absorbers at R± ~ 200 kpc, comparable 
to this R gmax (27). We assume the filamentary emission in the discovery spectrum extends 
from +50 kpc to +150 kpc at positive slit position, and —50 to —80 kpc at negative slit po¬ 
sitions, corresponding to an area of 15 sq. arcsec given our 1" slit (1130 kpc 2 ). The total area 
is equal to the Nq SO = 23 quasars observed to sufficient depth, times the 23 sq. arcsec area 
(1740 kpc 2 ) for which our slit intersects the annular region, implying a total area of 530 sq. 
arcsec (39922 kpc 2 ) surveyed. For reference, the total annular area around a single quasar is 
834 sq. arcsec (62832 kpc 2 ), so the effective area of our survey is 64% of the area around a 
single quasar. Hence we deduce a covering factor of filamentary emission of f c ,m = 0.028. 

Thus taken in aggregate, we find that only ~ 3% of the area between R± = 50 — 150 kpc 
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around quasars is covered by Lya emission with SB LyQ . > 1CT 17 erg s -1 cm -2 arcsec -2 , which 
can be compared to the much larger covering factor of 31% above this level, measured for 
emission over the same area from our NB image of SDSSJ0841+3921 (see Fig. la). The ag¬ 
gregate and individual covering factors can be commensurated if the majority of quasars do 
not exhibit large-scale Lya nebulae above SB Lya > 10 - 17 ergs -1 cm -2 arcsec -2 , but a small 
fraction ~ 10% have bright nebulae covering a large area, as in SDSSJ0841+3921, and the 
product of this ~ 10% frequency and ~ 31% covering factor set the aggregate ~ 3% effective 
covering factor for the entire population. This argument is consistent with our intuition that, 
given the random orientation of the slit, we could probably only have identified the asymmetric 
extended emission around SDSSJ0841+3921 for ~ 50% of the possible slit orientations, and 
given that 23 objects were surveyed to the required depth, this again suggests that ~ 10% of 
quasars exhibit giant nebulae. 

An independent estimate for the frequency of giant Lya nebulae around quasars can also be 
determined from a narrow band imaging survey conducted by our group. As part of a continuing 
program to search for fluorescent Lya emission powered by luminous z ~ 2 quasars, to date we 
have obtained deep narrow band imaging of a sample of 26 objects. These survey observations 
are analogous to the narrow-band imaging of SDSSJ 0841+3921 discussed in Supp. 1.2, in that 
they use custom designed narrow-band filters, however the targets constitute an independent 
sample of quasars, which are not part of the QPQ survey of (22). These observations have been 
conducted on multiple telescopes, which we briefly summarize. For the present purposes we 
define a giant Lya nebula to be extended emission on a scale > 50 kpc from the quasar with 
average SB Lya > 10 -17 ergs -1 cm -2 arcsec -2 . 

The quasar HE0109—3518 at z = 2.41 was observed with the Very Large Telescope/ -FOcal 
Reducer Spectrograph (VLT-FORS) using a custom designed narrow band filter (A cent er = 
4145A, FWHM a = 40A) for a total exposure time of about 20 hours, achieving a sensitiv¬ 
ity limit (la) of SBL ya = 8 x 10 -19 ergs -1 cm -2 arcsec -2 (61). No extended Lya nebula 
was discovered. We have imaged a total of eight quasars with the Keck LRIS imaging spec¬ 
trometer, of which three were observed with a filter centered on Lya redshifted to z — 2.279 
(Acenter = 3986A, FWHM A = 30A), and five with one centered at z = 2.190 (A cen ter = 3878A, 
FWHM a = 30A). The exposure times for the LRIS observations vary from 2 to 10 hours, re¬ 
sulting in sensitivity limits (la) ranging from SBL ya = 1.2 x 10 -18 ergs -1 cm -2 arcsec -2 to 
5 x 10 -19 . These VLT/LORS and Keck/LRIS observations specifically targeted hyper-luminous 
quasars V < 17, with the goal of understanding fluorescent Lyct emission from dark-galaxies 
(61) and the quasar CGM (41), and led to the discovery of a giant Lya nebula around the 
quasar UM287 at z — 2.28. Nebulae were not detected around any of the other quasars 
surveyed. Linally, using a custom filter centered on Lya at z — 2.253 (A cen ter = 3955A, 
FWHM a = 30A), we have observed a total of 17 quasars using the Gemini Multi Object 
Spectrograph (GMOS) on the telescope Gemini South. For the GMOS observations, a subset 
of 3 quasars were observed with longer integrations, typically 5 hours and achieving a depth 
of (la) SB Lyo , = 1.5 x 10 - 18 ergs -1 cm - 2 arcsec -2 ; whereas the other 14 quasars were ob¬ 
served in a fast survey mode with typical exposure times of about 2 hours and depth SB Lya = 
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3 — 4 x 10 - 18 ergs -1 cm - 2 arcsec -2 . Although the instruments used and depths achieved 
are heterogeneous, all of the aforementioned observations reached sufficient depth to detect 
Lya nebulae with SB Lya > 10 - 17 ergs -1 cm -2 arcsec -2 . Taking all of these observations to¬ 
gether, we have observed a total of 26 quasars and detected a single giant Lya nebula around 
UM287 (41), implying that 1/26 = 0.04 of quasars host giant Lya nebulae. Adopting the 
la confidence interval appropriate for the Poisson distribution in the small number regime, we 
deduce that the frequency of occurrence of giant nebulae is in the range 3 — 9%. 

Finally, an independent narrow band imaging survey of quasars at z ~ 2.7 is being car¬ 
ried out as part of the Keck Baryonic Structure Survey (KBSS), and the first results on the 
distribution of LAEs around quasars were recently published by Trainor et al. (118). Based 
on these narrow band imaging observations, (39) published the discovery and additional ob¬ 
servations of a giant Lya nebula around the quasar HS1549+19 with properties comparable to 
that around SDSSJ0841+3921. Trainor et al. quote a typical narrow-band imaging depth of 
m NB (3a) 26.7 for point sources, resulting from 5-7hr integrations on Keck LRIS using filters 
with FWHM r^j 80A. This corresponds to a (la) SBL ya — 1.5 x 10 18 ergs 1 cm 2 arcsec 2 , 
which although shallower than our own Keck LRIS survey (owing to their wider filters), it 
is nevertheless sufficient to detect giant nebulae with S 11 r jya > 10 -17 ergs -1 cm -2 arcsec -2 . 
There are 8 objects observed in the Trainor et al. sample, but it is unknown to us whether ad¬ 
ditional Lya nebulae are present, besides that around HS 1549+19. As such, one can place a 
lower limit on the frequency ofgiantLya nebulae to be 1/8 = 0.125, ora la Poisson confidence 
interval of 10 — 29%. 

To summarize, three independent surveys indicate that the frequency of occurrence of giant 
Lya nebulae around quasars is consistent with being ~ 10%, although there are still limitations 
due to relatively small samples, heterogeneity of the datasets, and the differences in the search 
techniques (narrow-band imaging versus spectroscopy). 


8 The Environments of Other Quasars with Giant Lyci Neb¬ 
ulae and the Impact of Lya Fluorescence 

Although quasar clustering indicates that the majority of z ~ 2 quasars reside in moderate 
overdensities, we speculate that ~ 10 % trace much more massive structures, and it is this subset 
which also exhibits giant Lya nebulae. SDSSJ0841+3921 clearly supports this interpretation, 
but there are indeed other examples. 

Of course, we also have NB imaging data for UM287, the quasar at z — 2.27 with a gi¬ 
ant Lya nebula recently discovered by our group (41). However, because UM287 is a hyper- 
luminous quasar, we cannot characterize its environment via the number counts of LAEs as we 
have for SDSSJ0841+3921. The reason for this is that distribution of LAEs on Mpc scales 
from hyper-luminous quasars is likely to be dramatically enhanced by Lya fluorescence (61), 
and it is not straightforward to disentangle the fluorescence and clustering enhancements. For 
example, consider the hyper-luminous quasar HE0109—3518, the proto-typical case for Lya 
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fluorescence by dark-galaxies studied by (61). The number of LAEs within R < Ah~ l cMpc 
of HE0109—3518 is ~ 8.5 times larger than the cosmic background number expected from the 
LAE luminosity function, whereas the quasar-LBG correlation function (a reasonable proxy for 
the quasar-LAE clustering) averaged over the volume probed by these narrow-band observa¬ 
tions would only predict an enhancement of ~ 1.5. One comes to similar conclusions about the 
number of LAEs discovered around the hyper-luminous quasars in the KBSS sample of (118). 
Progress on understanding the environment of UM287, and its relationship to the giant Lya 
nebulae, will require a sample of galaxies selected via broad-band photometry (i.e. LBG selec¬ 
tion), to avoid contamination from fluorescence. Nevertheless, while fluorescence precludes a 
simple interpretation of LAE clustering around UM287, we note that, intriguingly, follow-up 
observations also revealed a faint V ~ 22 companion quasar, indicating an overdensity of AGN 
similar to SDSSJ0841+3921. 

Note however that enhanced LAE clustering around SDSSJ0841+3921 due to Lyo; fluores¬ 
cence is expected to be negligible, because brightest quasar in the system, namely the f/g quasar 
with V = 19.8 is an order of magnitude fainter than the hyper-luminous quasars which exhibit 
this fluorescence effect. Radiative transfer simulations of Lya fluorescence around quasars (61) 
supports this conclusion. Furthermore, the four AGN in SDSSJ0841+3921 clearly constitute a 
dramatic overdensity of AGN, and there is no way for Lya fluorescence to enhance AGN clus¬ 
tering. Note that the other three AGN do not significantly enhance the ionizing luminosity and 
hence cannot boost the fluorescence signal. The fainter Type-1 AGN AGN2 and AGN3 have 
V ~ 23 (see Table S2), and thus contribute negligibly to the total ionizing photon production 
compared to the f/g quasar. Although we do not detect the continuum from the accretion disk 
of the Type-2 AGN AGN1, it is detected in the WISE bands W1 and W2 (see Table S2). At 
z ~ 2 these bands probe rest-frame ~ 1 //in. which are likely still dominated by emission from 
the reddest part of the accretion disk. The fact that AGN 1 is > 1 magnitude fainter than the f/g 
quasar indicates that its UV ionizing continuum (although obscured from our perespective) is 
very unlikely to be brighter than that of the f/g quasar. 

The KBSS quasar HS1549+19 at z — 2.85, which has been narrow band imaged by Trainor 
et al. (118), harbors a giant Lya nebula (39) with SB and properties comparable to that around 
SDSSJ0841+3921. Extensive imaging and spectroscopy of HS1549+19 has demonstrated 
that it is clearly a dramatic protocluster, with a Mpc-scale overdensity of Lyman break galaxies 
(LBGs) a factor of ~ 3.5 times larger than that around typical quasars (10, 38). Although this 
quasar is hyper-luminous, Lyct fluorescence is not an issue because the protocluster overdensity 
is based on the clustering of broad-band selected LBGs, and not LAEs. 

Finally, one may have concerns that the dramatic enhancements of LAEs around HzRGs 
and LABs, quantified by our clustering analysis in Supp. 4 and shown in Figs. 2 and S5, could 
be partly due to Lya fluorescence. This is a valid concern, because there is compelling evidence 
that the HzRGs are obscured quasars (119), and their ionizing photons, although obscured from 
our perspective, surely shine in other directions, and are likely to power their giant Lya nebu¬ 
lae. Similar considerations hold for LABs, although a direct association with obscured AGN has 
been more challenging to establish (18,120-123). However, we think that it is unlikely that flu- 
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orescence plays a significant role in the strong clustering of LAEs around HzRGs and LABs for 
two reasons. First, there is no evidence that the unobscured counterparts of HzRGs and LABs 
(i.e. the Type-1 quasar shining in other directions) are as bright as the hyper-luminous V < 17 
quasars, and hence sufficient to power LAE fluorescence on Mpc scales. Furthermore, even if 
the HzRGs and LABs were shining brightly in unobscured directions, the LAE enhancement is 
expected to be significantly smaller than for unobscured quasars, because given the cylindrical 
volume probed by NB imaging observations, the illuminated volume will be much larger for an 
unobscured quasar pointing toward the observer and hence along the major axis of the cylindri¬ 
cal survey volume, than for an obscured quasar which illuminates a volume randomly oriented 
relative to our line-of-sight direction (61). 

Thus, although based on only a handful of objects, existing observations support a picture 
where ~ 10% of radio-quiet quasars are surrounded by giant Lya nebulae, which acts as a 
signpost for the presence of a protocluster. A similar connection between AGN activity, giant 
Lya nebulae, and protocluster overdensities has been suggested by past work on HzRGs and 
LABs. Our clustering analysis in Supp. 4 (see Figs. 2 and S5) demonstrates that these overden¬ 
sities are indeed very large and well in excess of the environment of radio-quiet quasars. These 
protocluster overdensities are real, and are not the result of Lya fluorescence. 


9 Absorption Line Analysis 

Given the experimental design of our QPQ survey, we are able to explore characteristics of 
the gas surrounding the quadruple AGN in unprecedented detail. Fig. 4 shows strong Hi Lya 
absorption coincident with the redshift of the protocluster. We measure a rest-frame equivalent 
width W\ = 3.75±0.05 A, offset by a velocity Sv = cAA/A ~ +650 km s” 1 from the f/g quasar 
systemic redshift. The absence of very strong damping wings limits the HI column density to 
be IVhi < 10 19 7 cm~ 2 . We also present our favored model which has N m = 10 19 2± °' 3 cm -2 , 
suggesting optically thick gas at the Lyman limit. It is possible, in the presence of significant 
line-blending, that the Nm value could be much lower, but we consider this scenario very 
unlikely, since it would imply a gas metallicity exceeding the solar abundance (see below). 

Fig. S7 presents spectra of the b/g quasar centered at the wavelengths of a series of strong, 
commonly observed UV resonance-line transitions (including the several transitions shown in 
Fig. 4 of the Main text) in the rest-frame of the f/g quasar (based on its near-IR systemic red¬ 
shift). The gas exhibits especially strong low-ion absorption (Ol 1302, Cll 1334) and moder¬ 
ately strong intermediate ion absorption (Nil 1083, Sim 1206). The latter absorption occurs 
within the Lya forest and could be partially blended with coincident IGM absorption, but we 
believe such contamination is modest because the line-profiles track one another as well as the 
other low-ion transitions. Surprisingly, there is no statistically significant absorption at the high- 
ion transitions of Si +3 or C +3 . Even without performing a detailed photoionization analysis, we 
can already conclude that the gas is not highly ionized, which further suggests the gas is opti¬ 
cally thick to ionizing radiation. One draws a similar conclusion from the relative strengths of 
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OI 1302 to Sill 1304, where the latter is observed to have larger equivalent width in systems 
that are highly ionized (e.g., (124,125)), contrary to our measurements for this system. 

For each of the absorption lines, we have estimated rest-frame equivalent widths with sim¬ 
ple boxcar integration. Ionic column densities were estimated from the metal-line transitions 
assuming the linear curve-of-growth approximation and adopting lower limits to positive de¬ 
tections and upper limits for non-detections (Table S6). Because neutral oxygen 0° undergoes 
charge-exchange reactions with atomic hydrogen, the ratio N 0 o /N m is very nearly equal to the 
number ratio O/H over a wide range of physical conditions and ionizing radiation fields (25). 
Therefore, we estimate a conservative lower limit to the oxygen metallicity [O/H] > —1. To 
estimate abundances of the other ions and to constrain physical conditions of the gas, we must 
evaluate the ionization state with photoionization models (see the next section). 

Our previous studies have shown that cool gas with strong metal-line absorption from lower 
ionization states is common in the environment surrounding z ~ 2 quasars (26, 28). In Fig. S8 
we show the distribution of equivalent widths for Hi Lya and Cll 1334 from the set of QPQ 
pairs with separations less than 200 kpc. The absorption strength for SDSSJ0841+3921 lies 
toward the upper end of the distribution but is not extreme. 

In Table S6 we include estimates for the abundances of the observed elements assuming 
ionization corrections from the photoionization model discussed below and N m = 10 19 2 cm -2 . 
The observations require >1/10 solar metallicity and several measurements suggest solar (or 
even super-solar) enrichment. 


10 Joint Constraints on Protocluster Gas in Absorption and 
Emission 

Our preferred mechanism for powering the giant Lya nebula in SDSSJ0841+3921 is fluo¬ 
rescent emission from gas in the CGM of the protocluster powered by the radiation from the 
f/g quasar, with the gas residing in a population of cool dense clouds which are optically thin 
(IVhi < 10 17 ' 5 ) to ionizing radiation. Because our absorption line observations of the b/g quasar 
provide constraints on the gas properties in the protocluster, we adopt a novel approach of 
modeling the emission and absorption properties of this gas simultaneously, using the Cloudy 
photoionization code (126). First, we describe our model for the distribution of the cool clouds, 
and introduce simple analytical models for the fluorescent emission to build intuition about the 
scaling with gas properties. Then, we consider a sequence of photoionization models of the gas 
intercepted by the b/g sightline, which is however undetected in emission. Our absorption line 
measurements and photoionization modeling imply that the total cloud column density must be 
log 10 i+H = 20.4 ± 0.4. Under the reasonable assumption that the Lyo; emitting gas in the neb¬ 
ula has the same Am as that detected in absorption, we then show that reproducing the emission 
level requires very high gas volume densities tir — 2.0 cm -3 . 
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10.1 Cool Cloud Model and Emission Estimates 


We consider an idealized model whereby the gas in the protocluster halo resides in a popula¬ 
tion of cool T ~ 10 4 K spherical gas clouds, which have a single uniform hydrogen volume 
density ?z H and hydrogen column density iV H . The clouds are uniformly distributed throughout 
a spherical halo of radius R, and their aggregate covering factor, defined to be an average over 
the sphere, is f c . Four parameters (-R,71 h ,-/Vh, /c) completely describe the distribution of the 
gas (see ( 22 ) for further details), for example the total cool gas mass can be written 
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We now consider fluorescent Lya emission from this population of cool gas clouds powered 
by a bright central quasar. There are two regimes where simple expressions can be obtained for 
the resulting emission, namely when the gas clouds are either optically thin N m < 10 17 5 cm -2 
or optically thick N m > 10 l7 " 5 cm -2 to Lyman continuum photons. For the optically thin case, 
the gas clouds are highly ionized, and the average Lya SB from the halo will be 
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( 12 ) 

Note that because the gas is already highly-ionized, the SB is independent of the luminosity of 
the central ionizing source provided that it is bright enough to maintain the gas as optically thin. 

If the clouds are optically thick to ionizing radiation, they will no longer emit Lya propor¬ 
tional to the volume that they occupy because of self-shielding effects. Instead, a thin highly 
ionized skin will develop around each cloud, and nearly all recombinations and resulting Lya 
photons will originate in this skin. The cloud will then behave like a ‘mirror’, converting a frac¬ 
tion 7 / thi ck = 0.66 of the ionizing photons it intercepts into Lyct photons (127). In this regime, 
it can be shown (see (22)) that the fluorescent SB, averaged over the halos is 
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To obtain the second equality in eqn ( [13] ), we assume that the quasar spectral energy distribu¬ 
tion obeys the power-law form L u = L l/LL (i//r' LL ) Q?uv , blueward of i'll and adopt a slope of 
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ctuv = —1-7 consistent with the measurements of (128). The quasar ionizing luminosity is then 
parameterized by L Uhh , the specific luminosity at the Lyman edge. As we describe in detail be¬ 
low, we estimate that f/g quasar has log 10 = 30.5, which is the value adopted in eqn. ( |T3| ). 
In the optically thick regime, we see that the SB now depends on the luminosity of the quasar, 
as well as on the covering factor fc, but is independent of the gas properties (i.e. n^andNu). 

The smooth morphology of the emission in the Lya nebula implies a covering factor of 
fc 2d 0.5. Even if the emitting clouds have angular sizes much smaller than the PSF of our 
NB imaging, the morphology of the nebular emission would appear much more clumpy for 
lower covering factors. We have explicitly verified this by generating simulated images of 
mock Lya emission distributions with the required SB and a range of covering factors fc, 
and then convolving them with our NB-imaging PSF (129). In what follows we will always 
assume a covering factor of fc = 0.5 for both the gas we detect in emission, as well as that 
detected in absorption. This value is well motivated: a) for the emission it is based on the 
smooth morphology of the nebula; b) for the absorption it is approximately the covering factor 
of optically thick absorbers on R < 200 kpc scales in the quasar CGM (27). 

The average SB of the giant Lya nebula at R ~ 100 kpc is SB LyQ! = 1.3 x 10 -17 
ergs -1 cm -2 arcsec -2 , which we compute via the emission from the area within our 2 a SB 
contour which intersects the the annulus R = [75,125] kpc. This is several times lower than the 
expectation from eqn. ( [13] ) for the optically thick regime. Thus our simple analytical scaling 
relations already appear to favor a situation where the clouds are optically thin. We return to 
this comparison below using photoionization models. 

10.2 Details of Cloudy Modeling 

The expressions in eqns. < [I2| ) and ( fl3| ) provide reasonable estimates for the fluorescent SB, but 
they neglect other sources of Lya emission such as collisional excitation (i.e. cooling radiation), 
and pumped or ‘scattered’ Lya- photons arising from the diffuse continuum produced by the gas 
itsel. Note that we do not attempt to model scattering of Lya photons produced by central f/g 
quasar itself. Radiative transfer simulations of radiation from a hyper-luminous quasar through 
a simulated gas distribution (41), have shown that the contribution of scattered Lya line photons 
from the quasar do not contribute significantly to the Lya surface brightness of the nebula. 
Given that this calculation was for a hyper-luminous quasar a factor of ~ 10 times brighter than 
the f/g quasar, scattered line photons from the quasar are expected to contribute negligibly to 
the nebular emission., and they also do not treat the temperature dependence of the emission. 
To fully account for these effects, we construct photoionization models for a slab of gas being 
illuminated by a quasar using the Cloudy software package (126). Cloudy predicts the emergent 
spectrum from the gas, allowing us to compute its SB, which we then dilute by the covering 
factor fc = 0.5. 

We model the quasar spectral-energy distribution (SED) using a composite quasar spectrum 
which has been corrected for IGM absorption (128). This IGM corrected composite is important 
because it allows us to relate the g-band magnitude of the f/g quasar to the specific luminosity 
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at the Lyman limit L IJ] ] . For energies greater than one Rydberg, we assume a power law form 
L v = L UhL (u / u LL ) avv and adopt a slope of «uv = —1-7, consistent with the measurements 
of (128). We determine the normalization L Jy| l by integrating the Lusso et al. (128) composite 
spectrum against the SDSS filter curve, and choosing the amplitude to give the correct g-band 
magnitude of the f/g quasar (see Table S2), which gives a value of log 10 L VLL = 30.5. We 
extend this UV power law to an energy of 30 Rydberg, at which point a slightly different power 
law is chosen a = —1.65, such that we obtain the correct value for the specific luminosity at 
2 keV L u ( 2keV) implied by measurements of a ox , defined to be L v (2 keV)/Lj,(2500 A) = 
(^ 2 kev/^ o) Q ° x . We adopt the value a ox = —1-5 measured by (130) for SDSS quasars 
of comparable luminosities. An X-ray slope of a x = — 1, which is flat in vf v is adopted 
in the interval of 2-100kev, and above 100 keV, we adopt a hard X-ray slope of o H x = —2. 
For the rest-frame optical to mid-IR part of the SED, we splice together the composite spectra 
of (128), (84), and (131). These assumptions about the SED are essentially the standard ones 
used in photoionization modeling of AGN (e.g. (132)). We also include the contribution from 
the radiation field of the ambient UV background radiation field using the spectrum of (133), 
but it is completely sub-dominant relative to the f/g quasar at the distances that we consider. 

These Cloudy models require as an input the ionizing photon number flux <F(r) (see 
eqn. (fl~4|)) at the surface of the cloud. Because this varies with radius as r” 2 , in principle mod¬ 
eling the halo emission would require integrating over a distinct photoionization model at each 
radius. For simplicity, we restrict attention to just a few radii. For our model of the optically 
thick absorber, we set r = 250 kpc, which we believe is the likely location of the absorbing 
gas for several reasons. First, the nebula itself has an end-to-end size of ~ 500 kpc roughly 
centered on the f/g quasar, such that r = 250 kpc is a reasonable radius for the absorber given 
the measured impact parameter R± = 176 kpc of the b/g quasar sightline. Second, given this 
impact parameter, the characteristic velocity separation of the absorbing gas from systematic 
Umax = 600 kms -1 implied by the metal-line absorption (see Fig. S7), we can compute the con¬ 
ditional probability distribution P(< r\R±, u max ) that the distance between the quasar and the 
absorber is less than r, which is determined by the quasar-LLS correlation function which we 
have measured (24, 27). For the problem at hand, we find that P(< r\R±,v ma , x ) = 0.43 for 
r = 250 kpc, again indicating that this is a reasonable distance (see a similar discussion in § 4.1 
of (25)). According to eqn. 0. the SB arising from such optically thick gas scales with <f>, 
and hence will decrease as r~ 2 . For the gas producing the Lya emission, we model three radii 
r = [50,100, 250] kpc. We will argue that the emitting gas in the nebula is optically thin, and in 
this regime the SB is independent of $ (see eqn. (121). Thus considering a few radii to describe 
the emission is a reasonable approximation provided the gas remains optically thin at the radii 
in question. 

The Cloudy model of the cool cloud is completely specified by the value of <f>, the SED 
shape, the metallicity of the gas Z, the hydrogen volume density n H , and a ‘stopping criterion’, 
which sets the total column density of the cloud. Because we assume the f/g quasar and UVB 
are the only sources of ionizing photons, <I> and the SED shapes are fixed. For the metallicity, 
we adopt a value of log (Z/Z Q ) = —0.5 consistent with our lower-limit (log (Z/Z Q ) > —1) on 
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the metallicity of the optically thick absorber detected in the b/g quasar sightline (see Supp. 9). 
We have confirmed that the results discussed below are insensitive to the assumed gas metal¬ 
licity. We have also verified that our results are insensitive to the exact value of the UV slope 
a = —1.7, within the 1 a error measured by Lusso et al. a = —1.7 ± 0.6. Our models thus 
constitute a one-dimensional sequence parameterized by the volume density of the clouds nn. 
Photoionization models are known to be self-similar in the ionization parameter U = T /riuc, 
which is the ratio of the number density of ionizing photons to hydrogen atoms. Since $ has 
been fixed, our sequence of models is thus also a sequence in U, spanning the relevant range of 
ionization conditions. Our stopping criterion depends on what is known about the gas that we 
are attempting to model, which we elaborate on next. 

10.3 Photoionization Model of the Absorbing Gas 

Here we compare photoionization models to the properties of the optically thick absorber in 
the spectrum of the b/g quasar, lying at at an impact parameter of R± = 176 kpc from the f/g 
quasar. We follow the standard approach (25) and focus our analysis on multiple ionization 
states of individual elements to avoid uncertainties related to intrinsic abundances. Specifically, 
our column density measurements in Table S 6 yield the following limits on ionic ratios: 
C+/c 3 + > 0.8, Si + / Si 3+ > 1.3, and Si + /Si 2+ < 0.7. In addition, we measured the neutral 
column density to be Am = 10 19 ' 2 cm -2 , which we use as the stopping criterion of the Cloudy 
model. In other words, Cloudy will continue to add slabs of material until this neutral column 
is reached, guaranteeing that we satisfy this constraint. Finally, we do not detect Lyo emission 
from the vicinity of the background sightline, allowing us to set an upper limit of SB L y Q < 5.1 x 
10 -18 ergs -1 cm -2 arcsec -2 from the cool gas. This SB limit was computed via simulations 
whereby a fake circular source with radius of 3" was placed at the location of the b/g quasar. 
This exercise showed that only sources with SB LyQ > 3 x SB lcr would be clearly detectable, 
where SB lcr = 1.7x 10 -18 ergs -1 cm -2 arcsec -2 is our 1 a SB limit for 1.0 sq. arcsec apertures. 
Note that the individual clouds intercepted by our sightline may have extremely small sizes 
~ 10 pc making it impossible to angularly resolve emission from a single cloud. However, we 
have argued that these clouds have a large covering factor of / c = 0.5 such that we would be 
able to detect their aggregate emission near the location of the b/g sightline (see Fig. 1 in (22) 
for more details). 

Our sequence of photoionization models is compared to the observational constraints in 
Fig. S9. The three different ionic ratios considered paint a consistent picture for the ioniza¬ 
tion state, and we constrain the ionization parameter to be log 10 U = —2.7 ± 0.3, where we 
conservatively include an error of 0.3 dex to account for both statistical and systematic uncer¬ 
tainties in the photoionization modeling. Accordingly, the neutral fraction is constrained to be 
log 10 xhi = —1-2 ± 0.3, and the total hydrogen column log 10 A H = 20.4 ± 0.4, where the er¬ 
ror on Ah includes the quadrature sum of modeling uncertainty in the neutral fraction, and the 
uncertainty in our measurement of Ahi (see Supp. 9). Plugging this determination of the total 
hydrogen column density into eqn. ([IT]) we obtain a total cool gas mass of M c = 2.4 x 10 11 M 0 
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within a spherical halo R = 250 kpc, or the range 1.0 x 10 11 M 0 < M c < 6.5 x 10 11 M 0 im¬ 
plied by the la measurement and modeling error on A r tI . Note that given the self-similarity of 
the photoionization models with ionization parameter U, these results are essentially indepen¬ 
dent of the value of $ that we have assumed in our models, and hence the distance of the gas 
from the quasar. Lastly, we have calculated the ionization corrections to the ionic column den¬ 
sities measured for the absorbing gas and provide limits on the elemental abundances implied 
by each metal absorption line measurement in Table S 6 . 

According to the photoionization models, if indeed the absorbing gas lies at a distance 
r = 250 kpc and is illuminated by the f/g quasar, the Lya emission would be SB LyQ ~ 
3 — 4 x 10 -18 erg s -1 cm -2 arcsec -2 , which lies just below our detection limit SBL yQ! = 5.1 x 
10 -18 ergs -1 cm -2 arcsec -2 , and is hence consistent with our non-detection of extended emis¬ 
sion near the b/g quasar sightline. Indeed, because the gas is optically thick, and essentially 
absorbs all ionizing photons that impinge on it, the SB LyQ! is nearly independent of the volume 
density, and as expected from eqn. ( [13] ), depends primarily on $ and our assumed covering fac¬ 
tor fc = 0.5. The weak dependence on n H arises because, at the lowest ionization parameters 
and hence highest neutral fractions, as much as ~ 20% of the Lyo is produced by collisional 
excitation (green curve labeled collisions in Fig. S9). This is because collisional excitation of 
neutral hydrogen requires a large neutral fraction, and it is also exponentially sensitive to tem¬ 
perature. The low ri]j < 10 - 2 cm 3 and high ionization parameter log 10 U > — 1 models are 
highly ionized xhi < 10 -3 , and have higher temperatures T > 17, 000 K, both of which tend 
to suppress collisional excitation. As a result of the very weak dependence of SB LyQ on ri\\ in 
the optically thick regime, our non-detection of emission does not allow us to put interesting 
constraints on the volume density n H for the absorbing gas. 

Note that the fact that our photoionization model predicts Lyo- emission from the absorbing 
gas just below our detection limit results from our assumed distance of r = 250 kpc, since 
for optically thick gas SB LyQ oc </> cx r -2 (see eqn. 13). Had we assumed a distance equal 
to the impact parameter R± = 176 kpc the predicted emission would be a factor of two higher 
SB LyQ ~ 10 -17 ergs -1 cm -2 arcsec -2 , and hence detectable. As such it is informative to briefly 
discuss all possible scenarios that would result in a non-detection of Lyct from the absorbing 
gas. These are: 1) the gas is illuminated by the f/g quasar, but it lies at a distance r > 250 kpc 
such that its fluorescent Lya emisssion is below our detection threshold; 2) the gas is illuminated 
by the quasar and is at a distance comparable to the impact parameter, but the covering factor of 
optically thick gas is much lower than we have assumed fc = 0.5; 3) the absorbing gas is not 
illuminated by the f/g quasar but is instead shadowed, perhaps because from its vantage point the 
quasar is extincted by the same obscuring medium invoked in AGN unification models (67,68). 

Unfortuately, our current observations cannot distinguish between these three scenarios, 
however we favor scenario 1), that the gas lies at a distance ~ 250 kpc from the quasar, as 
shown in Fig. S9. While it is possible that the covering factor of optically thick gas is indeed 
lower than the f c = 0.5 that we have assumed, we consider this unlikely for two reasons. 
First, this would be in conflict with the statistical properties of optically thick absorption for the 
aggregate population of quasars (27). Although, SDSSJ0841+3921 is atypical given that it is 
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a protocluster and hosts a giant nebula, we already argued that the smooth morphology of the 
giant nebula also argues for a high covering factor f c > 0.5 of the emitting gas, to prevent the 
emission from appearing too clumpy. Under the reasonable assumption that the absorbing and 
emitting gas have the same origin, one expects their covering factors to be comparable. 

At first glance, the possibility that the absorbing gas is shadowed from the quasar radiation 
seems quite plausible. Indeed, we have previously argued in the QPQ series that the anisotropic 
clustering of LLSs around quasars implies that the optically thick gas in the quasar CGM is 
often shadowed (22-24). However, we robustly determine the ionization parameter of the ab¬ 
sorbing gas to be log 10 U = —2.7 ± 0.3. If obscuration shadows the gas from the quasar, then 
the radiation field would be expected to be the UV background, which given our ionization 
parameter, would then imply a gas density riu ~ 10 - 3 cm~ 3 . In the next section, we will see 
that the detection of bright Lya emission implies that the density of the emitting gas is three 
orders of magnitude higher n H ~ 1 cm -3 . As indicated in the middle panel of Fig. S9, if the 
gas is illuminated by the quasar at a distance r = 250 kpc, then gas density implied by our ion¬ 
ization parameter is then % — 0.7 cm -3 , which is then consistent with the value deduced from 
modeling the emission. Thus again invoking the very reasonable assumption that the absorbing 
and emitting gas have the same origin, the fact that we infer a comparable riu for the absorb¬ 
ing gas as required for the emitting gas, seems to favor a scenario where the absorbing gas is 
illuminated by the f/g quasar, but lies at a distance r > 250 kpc, such that the fluorescent Lya 
is below our threshold. Future deeper Lya observations could test this hypothesis if they can 
reach SE>Ly a — 10 - 18 ergs _ 1 cm _ 2 arcsec -2 levels, such that one could then detect emission 
from the absorbing gas out to a distance as far as r ~ 500 kpc. 

10.4 Photoionization Model of the Emitting Gas 

Here we consider photoionization models for the emitting gas under the assumption that it has 
a comparable column density jV H to the gas we detect in absorption. In modeling the emission 
we consider gas at r = 100 kpc which is the characteristic scale of the emission nebula. We 
assume that the total hydrogen column density does not vary strongly between these distances. 
This hypothesis is supported by a photoionization modeling study of a large sample of quasar 
CGM absorbers, which indicates that on average the total column density varies weakly (if 
at all) with impact parameter on such scales. We run the same Cloudy photoionization models 
as described previously, but with the modification that the stopping criterion is now chosen to 
be that the total hydrogen column A'n = 10 20 4 cm -2 , deduced in the previous section. 

In Fig. S10 we show the results of our sequence of photoionization models for the Lya 
emitting gas. The black curve illustrates the dependence of SB LyQ , on volume density n H for 
our fiducial choice of the radiation field intensity, where the gas lies at a distance of r = 100 kpc. 
The upper x-axis indicates the corresponding ionization parameter U at this distance. The red 
and blue curves are for stronger and weaker radiation fields, corresponding to r = 50 kpc 
and r = 250 kpc, respectively. For low values of «h the gas remains optically thin A r ni < 
10 17 ' 5 cm -2 , and the SB Lya exhibits a nearly linear dependence on the volume density n H and 
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is nearly independent of the ionizing radiation intensity, as predicted by eqn. 0- The gas 
becomes self-shielding and optically thick once N m > 10 1 ' 5 cm -2 , at which point the SB LyQ 
plateaus at a single value depending only on <1> (since we have fixed the covering factor fc = 
0.5) and hence on distance r, as predicted by eqn. ( fT3| ). 

These models clearly indicate that, given the column density of cool material Nu = 10 20 4 
cm -2 , the volume density must be % — 2.0 cm” 3 , in order to reproduce the observed average 
emission level at 100 kpc of SB Lya = 1.3 x 10” 17 ergs -1 cm” 2 arcsec” 2 (indicated by the 
horizontal dashed line in Fig. S10). Provided that the emitting gas is optically thin N m < 10 17 5 , 
this conclusion is essentially indepdenent of the radiation field and hence the assumed distance 
of the gas from the f/g quasar. In this regard, the relatively flat radial SBi, ya profile observed for 
the giant nebula is readily explained. It naturally arises because order of magnitude variations 
in the radiation field, corresponding to gas at distances ranging from r = 50 — 200 kpc, produce 
nearly the same level of Lyo emission, provided that IVh does not vary significantly across the 
nebula, as we have assumed. 

A corollary of our absorption-line constraint on Nn — 10 2 ° 4 cm” 2 and our emission con¬ 
straint on tt-h — 2.0 cm” 3 is that the implied absorption path length through the cool clouds is 
^cloud = N K /n K rs-/ 40 pc, implying that the emitting gas is distributed in extremely compact 
dense clouds. 
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Fig. SI: Two dimensional spectrum plotted as a x-map for the three slits in Fig. 3. The lower and 
middle rows of each image show Xsky (sky-subtracted only) and Xsky+PSF (sky and PSF subtracted), 
respectively as defined in § |T] (see ( 22 ) for additional details). The upper row shows the smoothed map 
Xsmth? (sky and PSF subtracted, then smoothed), which is helpful for identifying extended emission. A 
symmetric Gaussian kernel (same spatial and spectral widths in pixels) was adopted, with dispersion 
<7 sm th = 100 km s ' 1 (FWF1M = 235 km s -1 ), and ler spatial smoothing of 0.65" (FWFLM^l.5")- In the 
absence of extended emission, the distribution of pixel values in the the sky and PSF subtracted images 
should be Gaussian with unit variance. The lower and middle rows are displayed with a linear stretch 
ranging from — 6a to 6 cr indicated by the color-bar at middle right. The upper row has a stretch from 
—6a to 16<7 is adopted, and a 7 value chosen to decrease the range of color at the highest significance, as 
indicated by the color-bar at upper right. The horizontal dashed lines indicates the spectral traces for the 
f/g and b/g quasar, and the vertical dashed line indicates An = 0, corresponding to the systemic redshift 
of the f/g quasar. Note that the b/g quasar spectrum in Slit3 (right) has much lower S/N ratio than in 
Slit2 (center) because of much poorer seeing. 
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Fig. S2: Optical spectra of the four AGN at the same redshift. From top to bottom: f/g quasar, 
AGN1, AGN2, and AGN3. The gaps in the coverage of the spectra of AGN1 and AGN2 are due to the 
LRIS configuration chosen to cover both Lyo, C IV, and He II. The spectrum of AGN3 does not cover 
the Lya: line due to the limited blue sensitivity of DEIMOS. Note the different scale on the y-axis to 
better show the different emission lines, which are marked with vertical dotted lines and labeled. 
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Fig. S3: Comparison of the Line-ratios of AGN1 to other Type-2 AGN. The line-ratios of AGN1 
(red square) are compared to other Type-2 AGN compiled from the literature in the C IV /He II vs C I v/ 
C III] plane (left) as well as the C III] /He II vs C III] /C II] plane (right). The circles are individual mea¬ 
surements of HzRGs (black) from the compilation of (80), and narrow-line X-ray sources (cyan) and 
Seyfert-2s (blue) from the compilation of (81). Triangles indicate measurements from the composite 
spectra of HzRGs (orange) from (70), the composite Type-2 AGN spectra (purple) from (82), who split 
their population into two samples above and below Lya EW of 63A, and a composite spectrum of mid-IR 
selected Type-2 AGN (green) from (83). The stars represent measurements from composite spectra of 
Type-1 quasars, based on the analysis of (magenta) (84), (red-magenta) (85), and (blue-magenta) (86). 
The line-ratios for AGN1 lie in the locus of points spanned by other Type-2 objects in the literature. 
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Fig. S4: Confidence regions on correlation function parameters from maximum likeklihood clus¬ 
tering analysis. The 1 ,2, and 3a confidence regions are illustrated by the black, blue, and red contours, 
respectively. The maximum likelihood value is ro = 29.3 /i -1 Mpc and 7 = 1.5, which is indicated 
by the white cross. Due to the low signal-to-noise ratio of the clustering measurement, and the form 
adopted for the correlation function £(r) = (r/ro ) -7 (such that amplitude and slope are not indepen¬ 
dent parameters), a strong degeneracy exists between ro and 7 , such that neither is individually well 
determined. 
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Fig. S5: Binned correlation function of LAEs around HzRGs and LABs. The dimensionless corre¬ 
lation function 7 defined by eqn. (|9|) is indicated by the data points, with la error bars, determined by 
taking the 16th and 84th percentiles of the distribution of x resulting from our bootstrap procedure. The 
solid red line indicates the best-fit protocluster-LAE cross correlation function (ro = 29.3 h~ 1 Mpc and 
7 = 1.5) determined by our unbinned maximum likelihood procedure. The gray shaded region indicates 
the la error on our cross-correlation function fit, implied by the erorrs on the parameters ro and 7 . This 
region is determined by evaluating 7 for each sample of the boostrap error distributions of ro and 7 , and 
taking the 16th to 84th percentile confidence region of 7 . Our best fit correlation function (red line), and 
its quoted errors (gray shaded) clearly provides an acceptable representation of the binned data points 
and their respective errors. 
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Fig. S6: Near-IR spectrum of the f/g quasar SDSSJ084158.47+392121.0 (top) and of the b/g quasar 
BOSSJ084159.26+392140.0 (bottom) obtained with the NIRI spectrometer on the Gemini-North 
telescope. Each quasar shows the detection of strong [O III] emission and modest H/3 emission. The 
[O iii]A5008A doublet lines were centroided, giving observed wavelength centroids of A 0 bs = 1-5231 ± 
0.0001/mi and 1.60954=0.00006/rm for the f/g and b/g quasars respectively. This gives systemic redshifts 
Zf g = 2.0418 for the f/g quasar and 2.21370 for the b/g quasar. The estimated uncertainty on this 
redshift measurements is approximately 50 km s , taking into account both the error in our wavelength 
calibration and the estimated uncertainty when one uses [O III] to assess the systemic redshift of a quasar. 
This uncertainty on the measurement of the redshift is far lower than the width of our custom narrow- 
band filter (~ 2600km s -1 ). 
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Fig. S7: One-dimensional spectra of the b/g quasar centered at the wavelengths of a series of 
commonly observed UV-line transitions in the rest frame of the f/g quasar. The gas exhibits espe¬ 
cially strong low-ion absorption (OlA1302, SillA1304, CllA1334) and moderately strong intermediate 
ion absorption (NIIA1083, Si IIIA1206). The latter absorption occurs within the Lya forest and could be 
partially blended with coincident IGM absorption, but such contamination should be modest given that 
the line-profiles track one another and also the low-ion transitions. In contrast, there is no statistically 
significant absorption corresponding to the SilvA1393 and C IVA1548 transitions indicating the gas is 
not highly ionized. See Table S6 for a summary of the properties of these metal absorption lines. 
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Fig. S8: Equivalent widths of cool gas absorption in the QPQ dataset. Measurements of the equiv¬ 
alent widths of the Cll 1334 transition against the measurement for Hi Lyct. The measurements are 
drawn from our QPQ survey (27, 28 ) and the pair that led to the serendipitous discovery of the quadruple 
AGN system is marked separately. The gas probed by the background quasar (at an impact parameter 
of 176 kpc) shows strong low-ion absorption characterstic of the full population. This indicates a highly 
enriched and optically thick gas. 
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Fig. S9: Photoionization model of the optically thick absorber in b/g quasar spectrum. The gas 

is assumed to be at r = 250 kpc and exposed to the radiation field of the f/g quasar. The models are 
a one-dimensional sequence parameterized by the hydrogen volume density ng on the x-axis, which 
given that the radiation field is fixed, is also a proxy for the ionization parameter U shown on the upper 
x-axis. Upper panel: Variation of the total hydrogen column Am with nn/U, given that we require the 
model to reproduce our measured neutral column density Ami = 10 19 ' 2 . The y-axis on the right shows 
the corresponding neutral fraction. Our measured log 10 U (middle panel) in turn implies log 10 Am = 
20.4 ± 0.4, indicated by the vertical dotted line. Middle Panel: Variation of three different ionic ratios 
(Si + /Si 2+ in red; Si + /Si 3+ in green; C + /C 3+ in blue) with nn/U. Our measured lower/upper limits on 
these ratios from Table S6 are illustrated by the colored arrows. We obtain a consistent photoioinization 
solution implying log 10 U = — 2.7T0.3, indicated by the vertical dotted line. Lower Panel: Predicted SB 
of extended Lya emission from the absorbing gas as a function of nn/U. The red and green curves show 
the respective contributions from recombinations and collisional excitation, and the black curve shows 
the total Lya SB. The horizontal dashed line with arrows indicates our detection limit for an extended 
source near the b/g quasar sightline SBL yQ < 5.1 x 10 -18 ergs -1 cm -2 arcsec -2 . If the absorbing gas 
is at r = 250 kpc as we have assumed, our models are marginally consistent with a non-detection of 
extended Lya emission at that location. If the gas lies at larger distances r > 250 kpc, the predicted SB 
will decreaes as r 
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Fig. S10: Photoionization model of the Ly a emitting nebular gas. The models assume a total hy- 

, deduced from our photoionization modeling of the optically 


- 1020.4 


cm 


-2 


drogen column density N-r = 10 
thick absorber in the b/g sightline (see Fig. S9and § 10.3). The gas is assumed to be exposed to the 


radiation field of the f/g quasar at distances of r = 50, 100, and 250 kpc, indicated by the red, black, 
and blue curves, respectively. The models arc parameterized by the hydrogen volume density n h on 
the x-axis. The ionization parameter U corresponding to the gas at r = 100 kpc (black curves) is de¬ 
noted on the upper x-axis. Upper panel: Variation of the neutral hydrogen column IVhi with «h. The 
y-axis on the right shows the corresponding neutral fraction. Lower panel: Variation of the Lya SB 
with «h- For low ng the gas remains optically thin JVhi < 10 175 cm -2 , and the SB SBL yQ varies 
linearly with «h and is nearly independent of the ionizing radiation intensity. At higher nn the gas 
begins self-shielding A r ni > io 17 - 5 cm 2 , and the SB plateaus at a single value depending only on <f> 
and hence on distance r. The models indicate that a value of n h — 2.0 cm -3 , indicated by the ver¬ 
tical dotted line, is required in order to reproduce the observed average emission level at 100 kpc of 
SBL yQ , = 1.3 x 10 -17 ergs -1 cm -2 arcsec -2 , indicated by the horizontal dashed line in Fig. S10. This 
conclusion is essentially independent of the radiation field and hence the assumed distance of the gas 
from the f/g quasar. 
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Name 

RA 

DEC 

z 

o z 

(kms -1 ) 

e 

(") 

R L 
fkpc) 

f/g QSO 

08 41 58.47 

+39 21 21.0 

2.0412 

50 

— 

— 

AGN1 

08 41 58.24 

+39 21 29.1 

2.055 

400 

8.6 

71 

AGN2 

08 41 58.66 

+39 21 14.7 

2.058 

700 

6.6 

55 

AGN 3 

08 41 59.10 

+39 21 04.6 

2.050 

700 

18.0 

150 

b/g QSO 

08 41 59.25 

+39 21 40.0 

2.2138 

50 

21.1 

176 


Table SI: AGN properties. The coordinates, redshift z, redshift error a z , angular 9 and transverse separa¬ 
tion Rj_ from the f/g quasar arc listed for the four AGN in the protocluster, as well as the b/g quasar. 
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Line FWHM Flux W A 

(kms -1 ) (10~ 17 ergs -1 cm -2 ) (A) 

Lya. 1575 ± 36 76.5 ± 1.5 > 118.2 

NV.. — —0.4 ± 1.4 — 

CIV.. — 1.3 ±0.8 — 

Hell. 1154 ±276 2.4 ± 0.4 >15.8 

CIII]. 1434 ±171 2.6 ± 0.3 22.6 ±2.2 

CII]. — 1.0 ±0.2 11.2 ±2.6 


Table S3: Emision line measurements for the Type-2 quasar AGN1. Line fluxes were measured from the 
one-dimensional spectrum of the AGN, and errors determined from the la noise vector. For lines de¬ 
tected at S/N > 3, rest-frame equivalent widths W\ are computed. If the continuum lies above the la 
noise vector these are listed as values, otherwise we quote lower limits. The FWHM for emission lines 
for which our discovery spectra had sufficient S/N for a measurement (Lya, Hell, and Cm]) are also 
listed. 
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Tag 

RA 

DEC 

/Lya 

fx 

w Lya 

R± 




(10 -17 ergs -1 cm -2 ) 

(10 -19 ergs — 1 cm — 2 A — 1 ) 

(A) 

(kpc) 


f/g QSO 

08:41:58.470 

+39:21:20.97 

1708.76 

774.23 


72.61 

0 

AGN 2* 

08:41:58.685 

+39:21:15.00 

41.13 

36.57 


36.98 

50 

* 

08:41:57.982 

+39:21:25.88 

12.07 

2.04 


195.0 

58 


08:41:58.872 

+39:21:14.45 

4.14 

3.38 


40.28 

62 

AGN 1* 

08:41:58.238 

+39:21:28.87 

210.08 

62.35 


110.80 

65 


08:41:57.535 

+39:21:17.72 

1.33 

0.41 


106.14 

88 

* 

08:41:58.543 

+39:21:34.03 

21.45 

14.42 


48.91 

102 

* 

08:41:58.262 

+39:21:37.57 

14.07 

3.97 


116.50 

130 

AGN 3* 

08:41:59.107 

+39:21:04.67 

39.81 

37.76 


34.67 

139 

* 

08:41:58.286 

+39:21:39.20 

4.33 

0.95 


148.98 

143 


08:41:58.238 

+39:20:57.05 

3.42 

5.10 


22.06 

187 


08:41:59.035 

+39:21:46.27 

1.72 

1.09 


52.02 

203 


08:41:57.019 

+39:21:48.45 

1.13 

0.46 


80.45 

251 

Target 1* 

08:41:58.824 

+39:21:57.42 

17.22 

5.86 


96.66 

285 


08:41:59.762 

+39:21:55.52 

0.69 

0.14 


157.38 

293 


08:41:55.754 

+39:21:46.27 

0.86 

0.99 


28.70 

314 


08:41:54.958 

+39:21:14.99 

1.30 

0.49 


87.57 

321 


08:41:56.198 

+39:21:55.51 

0.78 

0.82 


31.32 

338 


08:41:53.926 

+39:21:33.75 

0.72 

0.52 


45.52 

422 


08:41:54.982 

+39:20:41.54 

0.85 

0.65 


43.33 

440 


08:41:59.038 

+39:20:20.88 

1.63 

2.77 


19.42 

470 


08:41:59.225 

+39:20:21.16 

0.85 

1.10 


25.30 

470 


08:41:56.762 

+39:20:19.79 

0.54 

0.89 


20.02 

500 


08:42:02.412 

+39:20:30.94 

0.93 

1.49 


20.44 

527 


08:41:57.794 

+39:20:11.63 

1.12 

1.56 


23.56 

543 

* 

08:41:52.003 

+39:21:14.16 

5.31 

3.44 


50.73 

586 


08:41:52.142 

+39:21:48.16 

1.21 

0.82 


48.76 

609 


08:41:55.778 

+39:20:03.75 

0.63 

0.63 


33.07 

648 


08:42:02.225 

+39:20:09.46 

0.96 

0.30 


105.43 

651 


08:42:02.412 

+39:20:06.74 

0.99 

0.98 


32.99 

678 

* 

08:41:52.802 

+39:20:21.41 

7.38 

3.22 


75.25 

690 


08:42:05.930 

+39:22:09.65 

0.82 

1.34 


20.01 

772 


08:41:58.027 

+39:23:01.05 

0.66 

0.80 


27.08 

779 


08:41:52.099 

+39:20:11.62 

1.86 

1.28 


47.73 

788 


08:41:51.722 

+39:20:16.51 

0.74 

0.19 


126.03 

789 


08:42:07.313 

+39:21:23.42 

0.56 

0.91 


20.07 

798 


08:41:49.610 

+39:21:22.58 

0.55 

0.85 


21.21 

800 


08:42:05.486 

+39:22:24.88 

0.59 

0.65 


29.80 

805 


08:41:51.209 

+39:20:17.33 

0.68 

0.30 


75.60 

821 


08:41:51.701 

+39:20:08.90 

0.62 

0.68 


30.16 

829 


08:41:59.340 

+39:23:07.04 

0.61 

0.80 


25.07 

829 


08:42:07.829 

+39:21:32.66 

0.92 

1.10 


27.59 

849 


08:41:51.701 

+39:20:01.01 

0.92 

1.33 


22.82 

872 


08:41:50.482 

+39:20:10.25 

0.67 

0.66 


33.42 

907 


08:41:51.067 

+39:20:02.09 

0.96 

0.11 


282.90 

907 


08:41:54.461 

+39:23:11.92 

0.91 

0.97 


30.68 

936 


08:41:49.982 

+39:22:30.30 

0.66 

0.85 


25.42 

937 


08:42:08.441 

+39:22:05.84 

1.71 

1.33 


42.09 

964 


08:42:07.826 

+39:20:09.72 

0.73 

0.69 


34.54 

1010 


08:42:10.267 

+39:21:13.89 

0.70 

0.60 


38.54 

1065 


08:41:51.482 

+39:23:15.72 

0.62 

0.92 


22.36 

1093 


08:42:10.759 

+39:21:24.49 

0.66 

0.25 


87.20 

1109 


08:42:10.010 

+39:22:10.45 

2.21 

0.27 


268.26 

1110 


08:42:12.590 

+39:21:39.16 

2.09 

1.83 


37.44 

1281 


08:42:13.382 

+39:20:05.34 

11.69 

4.20 


91.60 

1468 


08:42:11.702 

+39:23:13.80 

1.11 

0.73 


49.85 

1481 


08:42:12.840 

+39:19:12.32 

0.82 

0.62 


43.35 

1638 


08:42:03.629 

+39:17:19.50 

2.27 

1.16 


64.47 

1935 


08:41:52.363 

+39:17:15.40 

1.35 

1.44 


30.75 

1988 


08:42:03.677 

+39:17:02.64 

1.22 

0.65 


61.92 

2064 


08:42:12.223 

+39:16:41.66 

2.61 

2.83 


30.30 

2502 

Table S4: 

Properties of the candidate LAEs around the f/g quasar. Here / Lya is 

the Lya line flux in 

10" 17 

erg s 1 cm 2 , f\ is the flux density in 

the continuum in 

units of 10 19 erg s 1 cm 

- 2 A" 1 

estimated from the 


V-bantl, Wi, ya is the rest-frame Lya equivalent width, and R± the proper projected distance from the f/g quasar in 
kpc. LAEs are selected using the procedure described in § 4.1, namely they are required to have Hj, yo > 20A and 
fhya > 5.4 x 10 -18 erg s -1 . Objects which were used to measure the overdensity profile around the f/g quasar in 
Fig. 2 are denoted by a * in the column Tag, where the name of the source is also indicated. These sources had to 
satisfy a brighter flux limit, namely fhya > 4.0 x 10“ 11 erg s -1 cm -2 , and lie within 2.08' of the f/g quasar. 
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Name 

Type 

RA 

DEC 

z 

Filter FWHM (A) 

EW limit (A) 

Reference 

MRC 2048-272 

HzRG 

20 51 03.5 

-27 03 04.1 

2.06 

73 

20 

1 

MRC 1138-262 

HzRG 

11 40 48.2 

-26 29 09.5 

2.16 

65 

20 

1 

LABd05 

LAB 

14 34 11.0 

+33 17 32.6 

2.66 

201 

40 

2 

MRC 0052241 

HzRG 

00 54 29.8 

-23 51 31.1 

2.86 

66 

20 

1 

MRC 0943-242 

HzRG 

09 45 32.7 

-24 28 49.7 

2.92 

68 

20 

1 

LABI 

LAB 

22 17 26.0 

+00 12 37.7 

3.10 

80 

20 

3 

MRC 0316-257 

HzRG 

03 18 12.0 

-25 35 10.8 

3.13 

59 

20 

1 

TN J2009—3040 

HzRG 

20 09 48.1 

-30 40 07.4 

3.16 

59 

20 

1 


Table S5: Properties of objects analyzed in our protocluster-LAE correlation analysis. The nature of the 
source (HzRG or LAB), coordinates, redshift z, the FWHM of the NB filter used (A), the limiting EW 
of the LAE search (Wij m i t ), and the reference for the observations arc indicated. References used arc: 
(1) Venemans et al. 2007 (73); (2) Prescott et al. 2008 (79); (3) Nestor et al. 2011 (106). 
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Ion 

A 

(A) 

log / 

Vint 

(km s -1 ) 

Wa 

(mA) 

log TV 

[X/H] 

HI 

C 2 

1215.6701 

1334.5323 

-0.3805 

-0.8935 

[-225,1449] 

[355,1315] 

3750.0 ± 50.0 

1110.5 ±34.7 

19.2 ±0.3 

> 14.7 

> -1.19 

C 4 

1548.1949 

-0.7194 

[355,1315] 

< 331.5 

< 13.9 

< 0.72 


1550.7700 

-1.0213 

[355,1315] 

< 340.4 

< 14.2 

< 1.03 

N 2 

1083.9900 

-0.9867 

[480,850] 

216.3 ±31.8 

> 14.3 

> -1.00 

Ol 

1302.1685 

-1.3110 

[355,985] 

672.9 ± 14.0 

> 15.0 

> -0.95 

A12 

1670.7874 

0.2742 

[355,1315] 

335.0 ± 110.0 

> 12.9 

> -0.97 

Si 2 

1304.3702 

-1.0269 

[431,1315] 

259.3 ±17.5 

> 14.3 

> -0.81 

Si 3 

1206.5000 

0.2201 

[355,1315] 

951.0 ±43.6 

> 13.6 

> -0.80 

Si 4 

1393.7550 

-0.2774 

[355,1315] 

< 96.2 

< 13.0 

< -0.27 


1402.7700 

-0.5817 

[355,1315] 

< 93.7 

< 13.3 

< 0.02 


Table S6: Absorption line measurements. The columns indicate the Ion, rest wavelength A, oscillator 
strength log /, velocity interval n; nt used for the calculation of the rest equivalent widths, rest-equivalent 
width W\, ionic column density TV, and elemental abundances [X/H] implied by each metal absorption 
line measurement. The column densities TV reported assume the linear curve-of-growth approximation. 
The metal abundances [X/H] arc computed using ionization corrections from our favored Cloudy model 
presented in Section S10, and arc reported according to the standard convention of logarithmic relative 
to Solar (i.e. 0 = solar abundance). 
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